[0001]技术领域
[0002]本发明属于近海底质空间预测技术领域，具体涉及一种顾及不确定性的近海底质空间预测制图方法。
[0003]背景技术
[0004]近海底质类型图是近海海洋综合调查成果表达的基础图件之一，它在近海养殖、港口锚泊地选择、海洋牧场建设和近海沉积动力环境分析等方面都有着广泛的应用价值。高精度的近海底质类型图对于近海资源的科学开发利用至关重要。开展近海底质制图新方法研究，提升制图质量，既是海洋学科发展的需要，也是近海海洋资源管理与开发利用的需要，有着重要的现实意义。
[0005]底质类型是指根据沉积物中不同粒度组分(如黏土、粉砂、砂、砾石)的数量对比关系而对沉积物进行的一种分类和命名。目前国内外在底质制图领域普遍采用的底质分类和命名方案主要有Shepard方案和Folk方案两大类。Shepard方案一般只适用于无砾石沉积物的分类和命名，Folk分类方案可用于含砾和无砾沉积物的分类和命名，但它针对无砾石沉积物和含砾沉积物分别给出了两套分类和命名体系。
[0006]传统的底质类型图是海洋地质学家基于沉积物采样的粒度分析结果，结合水下地形、水动力条件和自身的专业知识人工勾绘而成。由于海洋地质学家的认知背景和知识经验存在差异，其勾绘过程总是带有各自的主观性，因而会在不同程度上对制图结果表达的客观性产生影响，并且人工勾绘效率低下。
[0007]近些年来，随着海底探测技术和地理信息技术的广泛应用，底质制图的定量化方法受到普遍重视。基于沉积物采样数据的空间预测制图方法是目前被广泛采用的定量化底质制图方法。它主要是通过对若干离散站位的底质沉积物进行观测或采样分析，获得沉积物的黏土、粉砂、砂、砾石等组分含量数据，借助空间统计和插值等方法来预测未采样站位处的各类底质组分的含量，然后再依据给定的底质分类和命名方案，对不同位置处的底质进行识别，进而实现区域底质类型的空间连续表达。该方法无需高密度采样，在大尺度范围的底质制图中有着非常显著的经费和时间成本优势；并且它对底质专题内容的表达可与现行的底质分类方案保持一致，其制图结果的沉积学内涵突出，并表现出较强的可比性，因而制图结果应用的普适性更为突出。
[0008]然而值得注意的是，当前基于沉积物粒度组分的底质空间预测制图方法有两个问题一直没有得到足够的重视：
[0009]一是沉积物粒度组分数据的成分数据属性对底质空间预测制图的影响。沉积物粒度组分数据属于单形空间(simplex space)的成分数据范畴，它具有“非负”(即各组分的取值均≥0)和“定和”(即各组分之和为定值)特性，这使得各组分之间存在闭合效应并由此产生伪相关，直接利用各组分数据进行空间预测时，难以保证各组分的空间预测结果总是能满足“定和”的要求，进而也无法开展底质类型识别。为解决这一问题，有学者提出运用对数比(Log-Ratio)转换方法来破解成分数据的闭合效应，进而实现成分数据由单形空间向欧氏空间(Euclidean space)的转换，并为经典统计学方法分析成分数据带来了方便。然而，对数比转换方法尽管解决了成分数据空间预测过程中的“定和”问题，但若成分数据中存在“0”值，则其应用会因“0”值不能取对数而受到限制。有学者认为沉积物粒度的空间分布往往是连续的，其组分的“0”值可能是由取样或分析测试结果的舍入误差引起的，因此建议用一个很小的非“0”值来替代成分数据中的“0”值，这样就可以解决成分数据中“0”值无法取对数的问题。但非“0”值替换也会带来一些问题，如在诸多的海洋底质调查中，无砾沉积物样品是普遍存在的，若对一些无砾沉积物样品人为地引入非“0”的极小值，会导致这些无砾沉积物的空间预测结果变成含砾沉积物，这在Folk分类方案中会出现原本应按照无砾沉积物进行分类和命名的位置转向按照含砾沉积物进行分类和命名。这显然对基于Folk分类方案的底质分类制图来说是不利的。
[0010]二是底质空间预测结果的不确定性对制图结果的影响。由于采样密度的空间差异和空间预测模型的不完备性是客观存在的，这使得空间预测模型对于不同空间位置处变量预测结果的方差是不一样的，也即预测结果的精度是不一样。而在目前底质空间预测制图方法中，不同位置底质类型的识别和命名主要是建立在该位置沉积物各粒度组分变量的预测结果之上的，而忽视了各粒度组分变量预测结果在该位置处的预测方差。这意味着，其底质制图的最终结果在不同空间位置处的可信度是不一样的，这对用图者来说是不利的。
[0011]发明内容
[0012]针对底质空间预测制图中存在的粒度组分数据的“定和”、“0”值和预测结果的不确定性问题，本发明提出一种顾及不确定性的底质空间预测制图方法，实现从概率测度的视角来开展底质类型的预测制图，并且它能充分顾及模型预测制图的不确定性。
[0013]为了实现上述目的，本发明采用的技术方案如下：
[0014]一种顾及不确定性的近海底质空间预测制图方法，其特征在于：方法如下：
[0015]S1：制图数据前处理，沉积物粒度组分数据的球坐标变换；
[0016]S2：空间预测，沉积物粒度组分转换变量的空间预测；
[0017]S3：多变量联合正态分布模型构建及底质类型识别，沉积物粒度组分转换变量空间预测结果的多变量联合正态模型构建及蒙特卡罗随机模拟和底质类型识别；
[0018]S4：制图结果表达，基于概率测度生成底质类型分布图及其不确定分布图。
[0019]进一步的，所述沉积物粒度组分实测数据是指黏土、粉砂、砂、砾石组分含量。
[0020]进一步的，S1中球坐标变换：沉积物粒度组分转换变量：(r,θ 2 ,…,θ p )，r为粒度组分之和的平方根；p为沉积物粒度组分数量。
[0021]进一步的，沉积物粒度组分转换变量预测结果μ i 、方差σ i 2 、不同变量预测结果间的协方差cov(θ i ,θ j )，i,j＝2,3,…,p。
[0022]进一步的，粒度组分转换变量预测结果的多变量联合正态分布表达：
[0023]N p-1 (μ,∑)
[0024]μ为各变量预测结果构成的向量；∑为各变量预测结果的方差和协方差构成的协方差阵。
[0025]进一步的，基于沉积物粒度组分转换变量预测结果的联合正态分布N p-1 (μ,∑)，运用蒙特卡罗随机模拟方法生成n组(如n＝5000)模拟值:
[0026]对每一组模拟值，利用球坐标变换的逆变换方法，将模拟结果转换为沉积物粒度组分在单形空间中的黏土、粉砂、砂和砾石组分的百分含量，据此再按照给定的底质分类方案，对该组模拟结果的底质类型进行识别；对于n组模拟结果，可采用同样的方法来对其底质类型进行识别，统计出不同类型底质出现占比，但n足够大时，该占比值可视为不同类型底质在指定位置处出现的概率大小；
[0027]运用上述方法遍历制图区域所有空间位置后，可获得给定底质分类方案中所有底质类型在整个制图区域出现的概率分布图。
[0028]进一步的，概率分布图硬化处理：通过比较制图区域内各空间位置处不同类型底质出现的概率大小，并按照概率最大准则将出现概率最高的那种底质类型赋予相应的空间位置，从而得到制图区域的底质类型分布图；
[0029]进一步的，概率分布图硬化处理过程会产生忽略不确定性，其大小可用信息熵或混淆指数来定量刻画，可以信息熵计算或混淆指数计算，具体如下：
[0030]信息熵的计算公式为：
[0031]式中H ij 为空间位置(i,j)处的不确定性信息熵，m为制图区域内所有底质类型的数量， 为空间位置(i,j)处第k类底质的出现的概率。H ij 的值在[0,1]之间，当H ij ＝0时，表明空间位置(i,j)处只出现某一种底质类型，硬化过程不会产生对其它底质类型的忽略，因此不会产生不确定性；当H ij ＝1时，表明在空间位置(i,j)处所有底质类型出现的概率是一样的，因此该位置归属任何底质类型都会带来最大的忽略不确定性；因此H ij 的值越大，不确定性就越高；但信息熵法只适用于概率值均不为0的情况；
[0032]混淆指数Confusion Index，CI的计算公式：
[0033]式中CI ij 为空间位置(i,j)处的混淆指数， 和 分别为空间位置(i,j)处不同类型底质出现概率值的最大值和第二大值；CI ij 的取值在[0,1]之间，其值越大，表明其不确定性越高。
[0034]通过上述技术方案可以得到以下有益效果：
[0035]制图数据的前处理提出了运用球坐标变换方法来破解沉积物粒度组分数据的闭合效应，它既满足的对沉积物粒度组分空间预测的需要(即满足预测结果的定和要求)，又不受“0”值约束，相较于传统的对数比转换方法，球坐标变换方法有着更为明显的普适性相较于传统的对数比转换方法，球坐标变换方法有着更为明显的普适性。
[0036]制图方法中提出根据粒度组分转换变量的空间预测结果及其方差和协方差来构建多变量联合正态分布模型，并用该分布模型来表达粒度组分转换变量空间预测结果的所有可能取值，相较于传统方法只重视变量预测结果，而忽视变量预测结果的不确定性(预测方差的差异性)和变量之间的关联性，新方法对粒度组分变量空间预测结果的表达更全面，它用统计分布来刻画预测结果，包含了预测结果本身所蕴含的不确定性。
[0037]对于制图区域底质类型图的生成是基于给定空间位置处出现概率最大的类别来确定的，也就是说，某种底质类型在给定空间位置处出现的概率最大，那么制图表达时，就将该类型底质赋予该位置。相对于传统底质制图方法，对某空间位置处底质类型的确定是根据该位置处沉积物粒度的黏土、粉砂、砂和砾组分含量来定义的，其结果是唯一的、确定性的。
[0038]总之：以球坐标转换方式来破解沉积物粒度组分数据的闭合效应，并且其转换过程不受“0”值约束，进而可实现对沉积物各粒度组分的有效预测(即满足“定和”要求)；以沉积物粒度组分转换变量的预测结果和方差为基础，建立粒度组分转换变量预测结果的联合正态分布模型，在此基础上运用蒙特卡罗方法随机生成若干组粒度组分变量预测结果的可能取值，经球坐标变换的逆变换后，可识别出每组模拟结果的底质类型，从而可统计出不同空间位置处不同底质类型发生的概率。由此实现从概率测度的视角来开展底质类型的预测制图，并且它能充分顾及模型预测制图的不确定性。
[0039]附图说明
[0040]图1是技术路线图。
[0041]具体实施方式
[0042]下面结合附图对本发明做进一步的说明：
[0043]如图1所示，本方法主要技术路线分四大阶段：(1)制图数据前处理——粒度组分数据的球坐标变换；(2)空间预测——粒度组分转换变量的空间预测；(3)多变量联合正态分布模型构建及底质类型识别——粒度组分转换变量空间预测结果的多变量联合正态模型构建及蒙特卡罗随机模拟和底质类型识别；(4)制图结果表达——制图区域底质类型专题图及其不确定性图的生成。
[0044]方法原理详述
[0045]1、沉积物粒度组分的球坐标变换
[0046]球坐标变换(Spherical Coordinate Transform，SCT)是一种将直角坐标系坐标转换为球面坐标系坐标的数据转换方法。假设某沉积物样品S是由p个粒度组分构成，即S＝(x 1 ,x 2 ,…,x ,p )，x i ≥0为粒度组分的百分含量，且 (i＝1,2,…,p)，将各组分含量分别做开平方的非线性变换，则有(x′ 1 ) 2 +(x′ 2 ) 2 +…+(x′ p ) 2 ＝10 2 ，其中x′ i ＝(x i ) 1/2 (i＝1,2,…,p)，那么S′＝(x′ 1 ,x′ 2 ,…,x′ p )∈R p 可被视为是半径r＝10的p维超球面上的一个点，数组(x′ 1 ,x′ 2 ,…,x′ p )即为该点在直角坐标系下的坐标值，相应地，该点也必然存在一个与之相对应的球坐标系下的坐标值。据此，我们可以将沉积物粒度组分数据S＝(x 1 ,x 2 ,…,x p )转换为p维球面坐标系下的球坐标(r,θ 2 ,…,θ p )∈Θ p ，其一般转换公式为：
[0047]与之相对应的逆变换公式为：
[0048]球坐标变换方法不仅可以破解沉积物粒度组分数据的闭合效应(即“定和”特性)，而且它不受成分数据中的“0”值约束，在成分数据的转换方面要较常见的对数比转换方法更具普适性。此外，沉积物粒度组分的球坐标变换结果在[0,π/2]之间，而对数比转换结果一般在(-∞，+∞)之间，这在一定程度上也降低了粒度组分转化结果(转换变量)的变异性，这对于建立更为稳健的空间统计分析模型是十分有利的。
[0049]2、沉积物粒度组分转换变量的空间预测
[0050]沉积物粒度组分变量属于区域化变量范畴，其空间分布往往具有一定的自相关性，因此可以借助地统计学方法来对其空间分布进行空间预测。由于沉积物粒度组分之间存在伪相关性，其各组分的转换变量之间通常也存在不同程度的相关性，因此本方法利用地统计学中的Cokriging方法来对沉积物粒度组分的转换变量进行协同估计。有关多变量Cokriging估计方法的原理很多文献都有详细的阐述，这里只对其作一简单介绍。
[0051]设沉积物粒度组分的K个转换变量构成的集合为{Z k (s)，k＝1,2,…,K}，s为空间位置的坐标。在满足二阶平稳假设条件下，交互半方差函数为
[0052]设Z p (s)为粒度组分转换变量{Z k (s)，k＝1,2,…,K}中的某一个待预测变量，也即主变量，其它为协变量，那么主变量Z p (s)在s 0 处的Cokriging预测值 为：
[0053]其中，n k 为预测点s 0 领域内的第k协变量采样点数据数量， 为分配给第k个协变量的第i个采样点的权重。根据最优无偏的假设条件，建立协同克里格的预测方程组：
[0054]式中，γ k′k (s i ，s j )是第k协变量和第k'协变量在位置s i 和s j 处的交互半方差值； 是位置s i 和待预测点s 0 处的交互半方差的平均值； 为第k变量的拉格朗日乘子。由上式可计算出权重λ，将其代入协同克里格的计算公式，可得到主变量 的预测值。
[0055]其相应的Cokriging预测方差为
[0056]3、粒度组分转换变量空间预测结果的蒙特卡罗模拟及其底质类型识别
[0057]对于任一空间位置，空间预测方法均可以获得该位置的变量预测结果(μ i )、预测方差(σ i 2 )和不同变量预测结果间的协方差cov(θ i ,θ j )，i,j＝2,3,…,p。因此，该位置沉积物粒度组分转换变量所有可能的预测结果可用(p-1)个转换变量的联合正态分布N p-1 (μ,∑)来刻画，其中，μ为各组分转换变量预测结果构成的向量，∑为各组分转换变量预测结果的方差和协方差构成的协方差阵。
[0058]基于沉积物粒度组分转换变量预测结果的联合正态分布N p-1 (μ,∑)，运用蒙特卡罗随机模拟方法生成n组(如n≥5000)模拟值:
[0059]对每一组模拟值，利用球坐标变换的逆变换方法，将模拟结果转换为沉积物粒度组分在单形空间中的黏土、粉砂、砂和砾石组分的百分含量，据此再按照给定的底质分类方案(如Folk方案、Shepard方案)对该组模拟结果的底质类型进行识别。对于n组模拟结果，可采用同样的方法来对其底质类型进行识别，在此基础上可进一步统计出不同类型底质出现占比，但n足够大时，该占比值可视为不同类型底质在指定位置处出现的概率大小。
[0060]运用上述方法遍历制图区域所有空间位置后，可获得给定底质分类方案中所有底质类型在整个制图区域出现的概率分布图。
[0061]4、底质类型图及其不确定性图的生成方法
[0062]底质类型图的生成可通过对制图区域内各类型底质概率分布图的“硬化”处理得到。“硬化”处理是指通过比较制图区域内各空间位置处不同类型底质出现的概率大小，并按照概率最大准则将出现概率最高的那种底质类型赋予相应的空间位置，从而得到制图区域的底质类型分布图。显然，该底质类型图是由制图区域所有空间位置处出现概率最高的底质类型构成。
[0063]由于底质类型图是基于概率最大准则得到了，因此当某一空间位置处的底质类型确定后，也必然忽略了其作为其它底质类型的可能性，因此“硬化”过程会产生忽略不确定性，其大小可用信息熵或混淆指数来定量刻画。
[0064]信息熵的计算公式为：
[0065]式中H ij 为空间位置(i,j)处的不确定性(信息熵)，m为制图区域内所有底质类型的数量， 为空间位置(i,j)处第k类底质的出现的概率。H ij 的值在[0,1]之间，当H ij ＝0时，表明空间位置(i,j)处只出现某一种底质类型，硬化过程不会产生对其它底质类型的忽略，因此不会产生不确定性；当H ij ＝1时，表明在空间位置(i,j)处所有底质类型出现的概率是一样的，因此该位置归属任何底质类型都会带来最大的忽略不确定性。因此H ij 的值越大，不确定性就越高。但信息熵法只适用于概率值均不为0的情况。
[0066]混淆指数(Confusion Index，CI)的计算公式：
[0067]式中CI ij 为空间位置(i,j)处的混淆指数， 和 分别为空间位置(i,j)处不同类型底质出现概率值的最大值和第二大值。CI ij 的取值在[0,1]之间，其值越大，表明其不确定性越高。
[0068]本发明具有如下三大创新点：
[0069]在制图数据的前处理方面
[0070]方法提出了运用球坐标变换方法来破解沉积物粒度组分数据的闭合效应，它既满足的对沉积物粒度组分空间预测的需要(即满足预测结果的定和要求)，又不受“0”值约束，相较于传统的对数比转换方法，球坐标变换方法有着更为明显的普适性。
[0071]在制图方法设计思想方面
[0072]方法从概率测度思想出发，提出根据粒度组分转换变量的空间预测结果及其方差和协方差来构建多变量联合正态分布模型，并用该分布模型来表达粒度组分转换变量空间预测结果的所有可能取值，相较于传统方法只重视变量预测结果，而忽视变量预测结果的不确定性(预测方差的差异性)和变量之间的关联性，新方法对粒度组分变量空间预测结果的表达更全面，它用统计分布来刻画预测结果，包含了预测结果本身所蕴含的不确定性。此外，新方法根据粒度组分转换变量预测结果的统计分布，运用蒙特卡罗方法随机模拟足量的可能取值，基于模拟值的球坐标变换的逆变换结果和给定底质分类方案开展底质类型识别，并统计不同底质类型的占比，也即不同类型在预测位置处发生的概率。因此，新方法对底质类型的预测结果不是基于粒度组分变量的唯一预测结果来确定的，而是根据粒度组分变量预测结果的所有可能取值，用概率测度的方法来确定，因此它顾及到了变量空间预测结果的不确定性。
[0073]在制图结果的表达方面
[0074]方法对于制图区域底质类型图的生成是基于给定空间位置处出现概率最大的类别来确定的，也就是说，某种底质类型在给定空间位置处出现的概率最大，那么制图表达时，就将该类型底质赋予该位置。相对于传统底质制图方法，对某空间位置处底质类型的确定是根据该位置处沉积物粒度的黏土、粉砂、砂和砾组分含量来定义的，其结果是唯一的、确定性的。因此，新方法从概率的视角来识别不同空间位置处的底质类型，其结果表达具有不确定性的内涵。此外，在制图结果的最终表达方面，传统制图方法只能提供一张制图区域的底质类型专题图，而新方法不但能提供底质类型的专题图，还能提供一张与专题图相匹配的不确定性大小分布图。
[0075]以上所述均为本发明的优选实施方式，对于本技术领域的普通技术人员，在不脱离本发明的原理前提下，对本发明的各种等价形式的修改均属于本申请所附权利要求的保护范围。
