CN103530904B - 一种基于克里金方法的水下地形数字高程建立方法 - Google Patents

Abstract

本发明公开了一种基于克里金方法的水下地形数字高程建立方法，其主要目的在于解决水深数据水深值变化幅度较大时引起数字高程建立过程中一些位置点存在水深误差较大的问题，所述方法步骤：初始水深数据的载入、根据水深数据建立水下不规则三角网地形、不规则三角网地形的空间优化、初始水深数据的变异结构分析、不规则三角网地形的水深克里金插值过程以及水深插值结果误差补偿。本发明方法通过不规则三角网地形的地形起伏选取克里金插值数据生成规则格网数字高程，相比直接通过初始水深数据插值生成规则格网数字高程，解决了水深数据水深值变化幅度较大引起规则格网数字高程误差较大的问题。

Description

一种基于克里金方法的水下地形数字高程建立方法

技术领域

本发明主要涉及水下地形数字高程模型领域，尤其涉及一种基于克里金方法的水下地形数字高程建立方法。

背景技术

现有的水下地形数字高程建立都是基于初始水深数据，直接通过插值方法完成的。水深地形数据库常直接应用于地形辅助导航，由于水下GPS信号不可达，因此导航的精确性要求水深地形数据库具备很高的精度。由于测深环境的复杂性使得水深数据存在系统误差和随机误差，进行粗差剔除后的初始水深数据呈不规则分布。直接从初始水深数据插值完成规则格网数字高程的过程中，存在插值方法、水深插值数据获取范围、水深插值结果修正的问题，其核心问题是待插值位置处水深数据范围的选取，无法根据地形变化特征动态选取水深数据，数据范围过小则不能正确的得出待插值位置水深值，数据范围过大，既有运算量大的问题，也会引入结果偏差。

因此，引入了一种基于克里金方法的水下地形数字高程建立方法，首先由初始水深数据完成不规则三角网地形，由于不规则三角网具有表达地形变化的特性，根据待插值位置附近不规则三角网地形的地形起伏选取合适步长的数据作为该位置处的水深插值数据，解决了数据选取的问题，基于克里金方法的线性加权水深值估计，采用权值修正和交叉验证的水深插值误差补偿，有益于提高规则格网数字高程模型的精度。

发明内容

本发明的目的在于：针对现有方法直接由初始水深数据通过插值方法完成水下地形数字高程的建立存在着水深插值数据选取不合理、插值结果没有误差补偿的问题，提出了一种基于克里金方法的水下地形数字高程建立方法，由本发明建立的水下地形数字高程具有精度高的优点。

本发明所述方法是通过完成初始水深数据的不规则三角网地形，根据地形起伏特征选取合适步长的水深数据进行普通克里金插值运算，再通过权值修正和交叉验证来完成水深插值结果的误差补偿，其技术方案如下：

一种基于克里金方法的水下地形数字高程建立方法，具体步骤如下：

步骤1采用公知的多波束水深测量方法测量水深值，获得水深数据

所述的水深数据表示经度为i、纬度为j的位置A_ij处的水深值，

步骤2建立不规则三角网

由初始水深数据在地球坐标系OENU中的NOE平面上的投影点相互连接建立不规则三角网，地球坐标系OE轴正向为地理东向，ON轴正向为地理北向，OU轴正向为天向：

步骤2.1确定投影点的凸壳

选取投影点中纬度值最小的投影点为起点，当存在多个投影点纬度值相等且最小时，取其中经度最小的投影点为起始顶点，记为X₀，由起始顶点X₀向其它投影点且A_ij≠0连线，取其中与X₀连线和OE轴夹角最小的投影点为第一顶点X₁，当存在多个这样的投影点时，取和起始顶点X₀距离值最大的投影点为第一顶点X₁，然后由第一顶点X₁向投影点且A_ij≠0，1连线，取其中与X₁连线和直线X₀X₁夹角最小的投影点为第二顶点X₂，以此类推，选取与X_k的连线和直线X_k-1X_k夹角最小的投影点为第k+1顶点X_k+1，k＝2，3，…n-2，直到所找到的第n顶点X_n和起始顶点X₀重合为止，n为所选取得到的顶点个数，由顶点X₀、X₁、X₂…X_n-1依次相连形成凸多边形，并以所述凸多边形为凸壳，顶点X₀、X₁、X₂…X_n-1为凸壳点，

步骤2.2由凸壳建立不规则三角网：

首先，由起始顶点X₀向凸壳的其它顶点依次连线，形成凸壳三角网，然后，对位于凸壳内的投影点作如下连接处理：

任选一个凸壳内的投影点，将所选投影点分别与包围所选投影点的三角形的顶点连线，遍历所有凸壳内的所有投影点，得到由三角形拼接而成的不规则三角网，

步骤2.3根据投影点与水深数据以及由三角形拼接而成的不规则三角网，构建由三角地形单元组成的不规则三角网地形，使得不规则三角网地形与不规则三角网形成对应关系，三角地形单元与不规则三角网中的三角形形成对应关系，步骤3不规则三角网地形空间优化

步骤3.1以任一共边的两个三角地形单元的所有顶点对应的水深数据点作为空间四面体的四个顶点构建空间四面体，所述空间四面体包括共边的两个三角地形单元及新增的两个三角地形单元，

步骤3.2如果空间四面体中的共边的两个三角地形单元的内角标准差之和小于新增的两个三角地形单元的内角标准差之和，则舍弃共边的两个三角地形单元并保留新增的两个三角地形单元；否则，舍弃新增的两个三角地形单元并保留共边的两个三角地形单元，

步骤3.3遍历所有共边的两个三角地形单元，得到空间优化的不规则三角网地形，

步骤4采用区域化变量分析方法，对初始水深数据进行结构变异分析，得出变异函数值表达式，所述的区域化变量分析方法为：

步骤4.1分别计算东西向、南北向、东北-西南向各距离的实际变异函数值，

在投影点中，以距离最小的一对投影点的距离值为一个距离单位，距离最大的一对投影点的距离值为h_max，分别筛选出投影点对之间的距离为h′_l＝l×h±Δh且过所述投影点对的直线与OE轴夹角θ′＝θ±Δθ的所有投影点对，计算筛选出的投影点对所在位置的水深数据对在计算距离为lh的变异函数值，h′_l为距离实际值，h为距离计算值且h＝0.5距离单位～1距离单位，l为距离因子，表示对向下取整运算，Δh为距离容限值且Δh＝0～0.05h，θ′为角度实际值，θ为角度计算值且东西向、南北向、东北-西南向分别取值为θ＝0°、θ＝90°、θ＝45°，Δθ为角度容限值且Δθ＝0°～2.5°，实际变异函数值计算如下：

${\mathrm{\γ}}_l=\frac{1}{2N_l}\underset{m=1}{\overset{N_l}{\mathrm{\Σ}}}{[Z_{A_{\mathrm{ij}}}-Z_{A_{\mathrm{ij}}^l}^m]}^2---\left(1\right)$

式中，γ_l为距离值为lh的水深数据对实际变异函数值，N_l为距离值为lh的水深数据对的个数，m为距离值为lh的水深数据对的序号，和分别为位置A_ij处及与其距离为lh的位置处的第m个水深数据对的一个水深值；

步骤4.2拟合实际变异函数值，得到三个方向的变异函数

所述拟合的方法为：选取球型理论变异函数，

${\mathrm{\&gamma;}}_{h^{\&prime;\&prime;}}=\left\{\begin{array}{cc}{C}_0,& h^{\&prime;\&prime;}=0\\ C_0+C(\frac{3}{2}\&CenterDot;\frac{h^{\&prime;\&prime;}}{a}-\frac{1}{2}\&CenterDot;\frac{h^{\&prime;\&prime;3}}{a^3}),& 0<h^{\&prime;\&prime;}\&le;a\\ C_0+C,& h^{\&prime;\&prime;}>a\end{array}\right.---\left(2\right)$

对于东西向、南北向、东北-西南向，计算处对应的理论变异函数值γ_h″，并根据理论变异函数值γ_h″和实际变异函数值γ_l，由最小二乘拟合法，得到三个方向变程值a、块金常数C₀、拱高C，带入球型理论变异函数，得到三个方向的计算变异函数；

步骤4.3三个方向的计算变异函数套合，求取变异函数值计算式，

所述的套合过程如下：

步骤4.3.1坐标系旋转变换

比较东西向、南北向、东北-西南向三个方向的变程值，将原坐标系OE轴和ON轴顺时针旋转β角形成O′EN系，0°≤β<180°，使OEN坐标系中变程值a最大的方向变为O′EN系横轴正向，则由OEN系到O′EN系的旋转矩阵Q的表达式为：

$Q=\left(\begin{array}{cc}\cos \mathrm{\&beta;}& \sin \mathrm{\&beta;}\\ -\sin \mathrm{\&beta;}& \cos \mathrm{\&beta;}\end{array}\right)---\left(3\right)$

步骤4.3.2坐标系压缩变换

将坐标系O′EN压缩变换成坐标系O″EN，由O′EN系到O″EN系的压缩矩阵P的表达式为：

$P=\left(\begin{array}{cc}1& 0\\ 0& \frac{R_1}{R_2}\end{array}\right)---\left(4\right)$

式中，P为压缩矩阵，R₁，R₂分别为东西向、南北向、东北-西南向三个方向中的最大变程值、最小变程值，

步骤4.3.3求变异函数值计算式

套合后的变异函数值计算式为：

$h^{\&prime;\&prime;\&prime;}=\sqrt{F^T{Q}^T{P}^T\mathrm{PQF}},F=\left(\begin{array}{c}{i}_0-i_1\\ j_0-j_1\end{array}\right)---\left(5\right)$

${\mathrm{\&gamma;}}_{h^{\&prime;\&prime;\&prime;}}=\left\{\begin{array}{cc}{C_0}^{\&prime;},& h^{\&prime;\&prime;\&prime;}=0\\ {C_0}^{\&prime;}+C^{\&prime;}(\frac{3}{2}\&CenterDot;\frac{h^{\&prime;\&prime;\&prime;}}{R_1}-\frac{1}{2}\&CenterDot;\frac{h^{\&prime;\&prime;\&prime;3}}{{R_1}^3}),& 0<h^{\&prime;\&prime;\&prime;}\&le;R_1\\ {C_0}^{\&prime;}+C^{\&prime;},& h^{\&prime;\&prime;\&prime;}>R_1\end{array}\right.---\left(6\right)$

式中，i_。，j₀，i₁，j₁分别表示在OEN系中的经纬度坐标值，F为由到的向量阵，F^T，P^T，Q^T分别为F，P，Q的转置，h′″为在O″EN系中与之间的距离，γ_h′″，为该距离值h′″的套合后的变异函数值，C₀′，C′分别为东西向、南北向、东北-西南向三个变异函数的块金常数、拱高的算术平均值，

步骤5将优化后的不规则三角网地形转换为规则格网地形，转换过程如下：

步骤5.1设定规则格网地形网格间距：

在初始水深数据在OENU系的投影点中，取距离最小的一对投影点的距离值的1／20～1／10为规则格网地形的网格间距，

步骤5.2采用普通克里金插值方法和水深插值误差修正方法估计水深值：

步骤5.2.1采用基于步长的方法，选取普通克里金插值方法的插值数据，所述的基于步长的方法如下：

对于步长的定义为：在不规则三角网中，包围待插值位置的三角形为处的0步长三角形，0步长三角形的顶点为处0步长位置点，与处0步长三角形共边相邻的三角形为处1步长三角形，1步长三角形的非共边顶点为处1步长位置点，以此类推，可得待插值位置处g步长三角形和位置点，和步长三角形相对应的为步长三角地形单元，步长位置点所在的水深数据为步长数据，设定初始粗糙度M₀＝10～30，G为步长阈值，表示步长的最大值，初值为1，计算待插值位置处步长阈值为G＝1内的三角地形单元的地形粗糙度M_G：

$M_G=\frac{\underset{g=0}{\overset{G}{\mathrm{\&Sigma;}}}\underset{u_g=1}{\overset{U_g}{\mathrm{\&Sigma;}}}{S}_{u_g}}{\underset{g=0}{\overset{G}{\mathrm{\&Sigma;}}}\underset{u_g=1}{\overset{U_g}{\mathrm{\&Sigma;}}}{S}_{u_g}^{\&prime;}}---\left(7\right)$

式中，u_g为g步长三角形的三角形编号，g＝0，1，…G，U_g为g步长三角形的三角形总数，为第u_g个g步长三角地形单元的空间平面面积，为第u_g个g步长三角形的平面面积，g步长三角形的平面面积为对应的g步长三角地形单元的空间平面面积在OENU系的投影面积，当M_G>M₀，令G＝G+1，进行新的步长阈值G的不规则三角网地形粗糙度M_G的迭代计算，直至迭代结果M_G≤M₀，并选取当前步长阈值G内所有的步长三角地形单元的顶点作为待插值位置处普通克里金插值方法的水深插值数据；

步骤5.2.2普通克里金权系数的计算，

由步骤5.2.1得待插值位置处步长阈值G内的总数为w个水深插值数据根据式(5)，式(6)，计算得待插值位置与水深插值数据在OEN系的投影位置在O″EN系中的距离值及变异函数值一对水深插值数据在OEN系的投影位置在O″EN系中的距离值及变异函数值e′＝1，2，…w且e′≠e，带入普通克里金方程组，

$\left[\begin{array}{ccccc}{\mathrm{\&gamma;}}_{h_1^1}& {\mathrm{\&gamma;}}_{h_1^2}& ...& {\mathrm{\&gamma;}}_{h_1^w}& 1\\ {\mathrm{\&gamma;}}_{h_2^1}& {\mathrm{\&gamma;}}_{h_2^2}& ...& {\mathrm{\&gamma;}}_{h_2^w}& 1\\ .& .& ...& .& .\\ .& .& {\mathrm{\&gamma;}}_{h_e^{e^{\&prime;}}}& .& .\\ .& .& ...& .& .\\ {\mathrm{\&gamma;}}_{h_w^1}& {\mathrm{\&gamma;}}_{h_w^2}& ...& {\mathrm{\&gamma;}}_{h_w^w}& 1\\ 1& 1& ...& 1& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&lambda;}}_1\\ {\mathrm{\&lambda;}}_2\\ .\\ .\\ .\\ {\mathrm{\&lambda;}}_w\\ \mathrm{\&mu;}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&gamma;}}_{h_0^1}\\ {\mathrm{\&gamma;}}_{h_0^2}\\ .\\ .\\ .\\ {\mathrm{\&gamma;}}_{h_0^w}\\ 1\end{array}\right]---\left(8\right)$

由普通克里金方程组得克里金权系数λ_e，e＝1，2，…w；

步骤5.2.3对普通克里金权系数进行修正，

所述修正的方法为：对权值为负的λ_e置0，其余权值修正为：

${\mathrm{\&lambda;}}_e^{\&prime;}=\frac{{\mathrm{\&lambda;}}_e}{\underset{k=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_k}---\left(9\right)$

式中，λ_e为修正前的权值，λ_k为对负权值进行置0处理后的权值，λ′_e为修正后的权值，

步骤5.2.4普通克里金插值方法的水深插值结果及误差修正，

待插值位置处的水深插值及插值方差的计算：

${\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^0}=\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_e^{\&prime;}{Z}_{A_{\mathrm{ij}}^e}---\left(10\right)$

${\mathrm{\&sigma;}}_{A_{\mathrm{ij}}^0}=\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_e^{\&prime;}{\mathrm{\&gamma;}}_{h_0^e}+\mathrm{\&mu;}---\left(11\right)$

对待插值位置处每一个水深插值数据e＝1，2，…w进行水深值交叉验证，取水深插值数据e，＝1，2，…w，e′≠e，由式(5)、式(6)、式(7)、式(8)、式(9)、式(10)、式(11)计算得位置处的水深计算值及误差

修正后位置处的普通克里金插值水深值为：

${\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^0}^{\&prime;}={\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^0}+\frac{{\mathrm{\&sigma;}}_{A_{\mathrm{ij}}^0}}{\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&sigma;}}_{A_{\mathrm{ij}}^e}}\&CenterDot;\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}(Z_{A_{\mathrm{ij}}^e}-{\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^e})---\left(12\right)$

步骤6输出规则格网数字地形。

与现有技术相比，本发明具有如下优点：

1)通过初始水深数据构造不规则三角网地形，由不规则三角网地形的地形粗糙度和预设的地形粗糙度进行比较，当大于预设粗糙度时，表明地形变化剧烈，则增加数据选取范围直至选取数据的不规则三角网地形粗糙度不大于预设值，相比于现有插值方法中固定数据个数和固定数据范围的数据选取方式或者克里金方法中常用的基于变程的数据选取方式，本方案中克里金方法基于粗糙度和己知步长阈值内不规则三角地形单元的关系来动态的与选取相应步长阈值内的水深数据，既避免了数据选取范围过小引起克里金插值结果误差大，也避免了数据范围选取过大从而插值过程出现克里金方程组不可解以至于无法进行插值，提高构建规则格网地形模型的效率和精度；

2)通过对克里金方程组解得的负权值进行置零，从而消除了负权值的水深数据值对插值结果水深值产生的误差；

3)通过对选取的水深插值数据进行交叉验证，可得实际水深插值数据的计算值和真值的偏差，用这些水深插值数据的偏差可以准确衡量这些水深数据局部范围内普通克里金插值水深估计的误差，进而估计待插值位置水深值的计算偏差，从而误差修正后的水深差值结果精度更高。

附图说明

图1为基于克里金方法的水下数字高程模型建立方法流程图。

图2凸壳三角网示意图。

图3为不规则三角网地形空间内角标准差优化一次优化过程原理图。

图4为初始水深数据结构变异分析流程图。

图5为坐标系旋转变换图。

图6为不规则三角网数字高程转化为规则格网数字高程流程图

图7为优化后的不规则三角网地形图

图8为基于步长的克里金插值数据选取示意图

图9为输出的规则格网地形图

具体实施方式

下面结合具体实施例和附图，进一步阐明本发明的技术方案，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域普通技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

一种基于克里金方法的水下地形数字高程建立方法，其特征为：

步骤1采用公知的多波束水深测量方法测量水深值，获得水深数据

所述的水深数据表示经度为i、纬度为j的位置A_ij处的水深值，步骤2建立不规则三角网

由初始水深数据在地球坐标系OENU中的NOE平面上的投影点相互连接建立不规则三角网，地球坐标系OE轴正向为地理东向，ON轴正向为地理北向，OU轴正向为天向：

步骤2.1确定投影点的凸壳

选取投影点中纬度值最小的投影点为起点，当存在多个投影点纬度值相等且最小时，取其中经度最小的投影点为起始顶点，记为X₀，由起始顶点X₀向其它投影点且A_ij≠0连线，取其中与X₀连线和OE轴夹角最小的投影点为第一顶点X₁，当存在多个这样的投影点时，取和起始顶点X₀距离值最大的投影点为第一顶点X₁，然后由第一顶点X₁向投影点且A_ij≠0，1连线，取其中与X₁连线和直线X₀X₁夹角最小的投影点为第二顶点X₂，以此类推，选取与X_k的连线和直线X_k-1X_k夹角最小的投影点为第k+1顶点X_k+1，k＝2，3，…n-2，直到所找到的第n顶点X_n和起始顶点X₀重合为止，n为所选取得到的顶点个数，由顶点X₀、X₁、X₂…X_n-1依次相连形成凸多边形，并以所述凸多边形为凸壳，顶点X₀、X₁、X₂…X_n-1为凸壳点，

步骤2.2由凸壳建立不规则三角网：

首先，由起始顶点X₀向凸壳的其它顶点依次连线，形成凸壳三角网，然后，对位于凸壳内的投影点作如下连接处理：

任选一个凸壳内的投影点，将所选投影点分别与包围所选投影点的三角形的顶点连线，遍历所有凸壳内的所有投影点，得到由三角形拼接而成的不规则三角网，

步骤2.3根据投影点与水深数据以及由三角形拼接而成的不规则三角网，构建由三角地形单元组成的不规则三角网地形，使得不规则三角网地形与不规则三角网形成对应关系，三角地形单元与不规则三角网中的三角形形成对应关系，步骤3不规则三角网地形空间优化

步骤3.1以任一共边的两个三角地形单元的所有顶点对应的水深数据点作为空间四面体的四个顶点构建空间四面体，所述空间四面体包括共边的两个三角地形单元及新增的两个三角地形单元，

步骤3.2如果空间四面体中的共边的两个三角地形单元的内角标准差之和小于新增的两个三角地形单元的内角标准差之和，则舍弃共边的两个三角地形单元并保留新增的两个三角地形单元；否则，舍弃新增的两个三角地形单元并保留共边的两个三角地形单元，

步骤3.3遍历所有共边的两个三角地形单元，得到空间优化的不规则三角网地形，

步骤4采用区域化变量分析方法，对初始水深数据进行结构变异分析，得出变异函数值表达式，所述的区域化变量分析方法为：

步骤4.1分别计算东西向、南北向、东北-西南向各距离的实际变异函数值，

在投影点中，以距离最小的一对投影点的距离值为一个距离单位，距离最大的一对投影点的距离值为h_max，分别筛选出投影点对之间的距离为h′_l＝l×h±Δh且过所述投影点对的直线与OE轴夹角θ′＝θ±Δθ的所有投影点对，计算筛选出的投影点对所在位置的水深数据对在计算距离为lh的变异函数值，h′_l为距离实际值，h为距离计算值且h＝0.5距离单位～1距离单位，l为距离因子，表示对向下取整运算，△h为距离容限值且Δh＝0～0.05h，θ′为角度实际值，θ为角度计算值且东西向、南北向、东北-西南向分别取值为θ＝0°、θ＝90°、θ＝45。，Δθ为角度容限值且Δθ＝0°～2.5°，实际变异函数值计算如下：

${\mathrm{\&gamma;}}_l=\frac{1}{{2N}_l}\underset{m=1}{\overset{N_l}{\mathrm{\&Sigma;}}}{[Z_{A_{\mathrm{ij}}}-Z_{A_{\mathrm{ij}}^l}^m]}^2---\left(1\right)$

式中，γ_l为距离值为lh的水深数据对实际变异函数值，N_l为距离值为lh的水深数据对的个数，m为距离值为lh的水深数据对的序号，和分别为位置A_ij处及与其距离为lh的位置处的第m个水深数据对的一个水深值；

步骤4.2拟合实际变异函数值，得到三个方向的变异函数

所述拟合的方法为：选取球型理论变异函数，

${\mathrm{\&gamma;}}_{h^{\&prime;\&prime;}}=\left\{\begin{array}{cc}{C}_0,& h^{\&prime;\&prime;}=0\\ C_0+C(\frac{3}{2}\&CenterDot;\frac{h^{\&prime;\&prime;}}{a}-\frac{1}{2}\&CenterDot;\frac{h^{\&prime;\&prime;3}}{a^3}),& 0<h^{\&prime;\&prime;}\&le;a\\ C_0+C,& h^{\&prime;\&prime;}>a\end{array}\right.---\left(2\right)$

对于东西向、南北向、东北-西南向，计算处对应的理论变异函数值γ_h″，并根据理论变异函数值γ_h″和实际变异函数值γ_l，由最小二乘拟合法，得到三个方向变程值a、块金常数C₀、拱高C，带入球型理论变异函数，得到三个方向的计算变异函数；

步骤4.3三个方向的计算变异函数套合，求取变异函数值计算式，

所述的套合过程如下：

步骤4.3.1坐标系旋转变换

比较东西向、南北向、东北-西南向三个方向的变程值，将原坐标系OE轴和ON轴顺时针旋转β角形成O′EN系，0°≤β<180°，使OEN坐标系中变程值a最大的方向变为O′EN系横轴正向，则由OEN系到O′EN系的旋转矩阵Q的表达式为：

$Q=\left(\begin{array}{cc}\cos \mathrm{\&beta;}& \sin \mathrm{\&beta;}\\ -\sin \mathrm{\&beta;}& \cos \mathrm{\&beta;}\end{array}\right)---\left(3\right)$

步骤4.3.2坐标系压缩变换

将坐标系O′EN压缩变换成坐标系O″EN，由O′EN系到O″EN系的压缩矩阵P的表达式为：

$P=\left(\begin{array}{cc}1& 0\\ 0& \frac{R_1}{R_2}\end{array}\right)---\left(4\right)$

式中，P为压缩矩阵，R₁，R₂分别为东西向、南北向、东北-西南向三个方向中的最大变程值、最小变程值，

步骤4.3.3求变异函数值计算式

套合后的变异函数值计算式为：

$h^{\&prime;\&prime;\&prime;}=\sqrt{F^T{Q}^T{P}^T\mathrm{PQF}},F=\left(\begin{array}{c}{i}_0-i_1\\ j_0-j_1\end{array}\right)---\left(5\right)$

${\mathrm{\&gamma;}}_{h^{\&prime;\&prime;\&prime;}}=\left\{\begin{array}{cc}{C_0}^{\&prime;},& h^{\&prime;\&prime;\&prime;}=0\\ {C_0}^{\&prime;}+C^{\&prime;}(\frac{3}{2}\&CenterDot;\frac{h^{\&prime;\&prime;\&prime;}}{R_1}-\frac{1}{2}\&CenterDot;\frac{h^{\&prime;\&prime;\&prime;3}}{{R_1}^3}),& 0<h^{\&prime;\&prime;\&prime;}\&le;R_1\\ {C_0}^{\&prime;}+C^{\&prime;},& h^{\&prime;\&prime;\&prime;}>R_1\end{array}\right.---\left(6\right)$

式中，i₀，j₀，i₁，j₁分别表示在OEN系中的经纬度坐标值，F为由到的向量阵，F^T，P^T，Q^T分别为F，P，Q的转置，h′″为在O″EN系中与之间的距离，γ_h′″为该距离值h′″的套合后的变异函数值，C₀′，C′分别为东西向、南北向、东北-西南向三个变异函数的块金常数、拱高的算术平均值，

步骤5将优化后的不规则三角网地形转换为规则格网地形，转换过程如下：

步骤5.1设定规则格网地形网格间距：

在初始水深数据在OENU系的投影点中，取距离最小的一对投影点的距离值的1／20～1／10为规则格网地形的网格间距，

步骤5.2采用普通克里金插值方法和水深插值误差修正方法估计水深值：

步骤5.2.1采用基于步长的方法，选取普通克里金插值方法的插值数据，所述的基于步长的方法如下：

对于步长的定义为：在不规则三角网中，包围待插值位置的三角形为处的0步长三角形，0步长三角形的顶点为处0步长位置点，与处0步长三角形共边相邻的三角形为处1步长三角形，1步长三角形的非共边顶点为处1步长位置点，以此类推，可得待插值位置处g步长三角形和位置点，和步长三角形相对应的为步长三角地形单元，步长位置点所在的水深数据为步长数据，设定初始粗糙度M₀＝10～30，G为步长阈值，表示步长的最大值，初值为1，计算待插值位置处步长阈值为G＝1内的三角地形单元的地形粗糙度M_G：

$M_G=\frac{\underset{g=0}{\overset{G}{\mathrm{\&Sigma;}}}\underset{u_g=1}{\overset{U_g}{\mathrm{\&Sigma;}}}{S}_{u_g}}{\underset{g=0}{\overset{G}{\mathrm{\&Sigma;}}}\underset{u_g=1}{\overset{U_g}{\mathrm{\&Sigma;}}}{S}_{u_g}^{\&prime;}}---\left(7\right)$

式中，u_g为g步长三角形的三角形编号，g＝0，1，…G，U_g为g步长三角形的三角形总数，为第u_g个g步长三角地形单元的空间平面面积，为第u_g个g步长三角形的平面面积，g步长三角形的平面面积为对应的g步长三角地形单元的空间平面面积在OENU系的投影面积，当M_G>M₀，令G＝G+1，进行新的步长阈值G的不规则三角网地形粗糙度M_G的迭代计算，直至迭代结果M_G≤M₀，并选取当前步长阈值G内所有的步长三角地形单元的顶点作为待插值位置处普通克里金插值方法的水深插值数据；

步骤5.2.2普通克里金权系数的计算，

由步骤5.2.1得待插值位置处步长阈值G内的总数为w个水深插值数据根据式(5)，式(6)，计算得待插值位置与水深插值数据在OEN系的投影位置在O″EN系中的距离值及变异函数值一对水深插值数据 在OEN系的投影位置 在O″EN系中的距离值及变异函数值e′＝1，2，…w且e′≠e，带入普通克里金方程组，

由普通克里金方程组得克里金权系数λ_e，e＝1，2，…w；

步骤5.2.3对普通克里金权系数进行修正，

所述修正的方法为：对权值为负的λ_e置0，其余权值修正为：

${\mathrm{\&lambda;}}_e^{\&prime;}=\frac{{\mathrm{\&lambda;}}_e}{\underset{k=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_k}---\left(9\right)$

式中，λ_e为修正前的权值，λ_k为对负权值进行置0处理后的权值，λ′_e为修正后的权值，

步骤5.2.4普通克里金插值方法的水深插值结果及误差修正，待插值位置处的水深插值及插值方差的计算：

${\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^0}=\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_e^{\&prime;}{Z}_{A_{\mathrm{ij}}^e}---\left(10\right)$

${\mathrm{\&sigma;}}_{A_{\mathrm{ij}}^0}=\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_e^{\&prime;}{\mathrm{\&gamma;}}_{h_0^e}+\mathrm{\&mu;}---\left(11\right)$

对待插值位置处每一个水深插值数据进行水深值交叉验证，取水深插值数据由式(5)、式(6)、式(7)、式(8)、式(9)、式(10)、式(11)计算得位置处的水深计算值及误差

修正后位置处的普通克里金插值水深值为：

${\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^0}^{\&prime;}={\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^0}+\frac{{\mathrm{\&sigma;}}_{A_{\mathrm{ij}}^0}}{\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}{\mathrm{\&sigma;}}_{A_{\mathrm{ij}}^e}}\&CenterDot;\underset{e=1}{\overset{w}{\mathrm{\&Sigma;}}}(Z_{A_{\mathrm{ij}}^e}-{\stackrel{\&OverBar;}{Z}}_{A_{\mathrm{ij}}^e})---\left(12\right)$

步骤6输出规则格网数字地形。

对本发明的有益效果进一步说明如下：

(1)水深插值数据的高效选取

仿真采用经度117.067-117.883，纬度33-34.35的1300个，最小间距为1850m的非等间距东海水深数据。由步骤4进行东西方向、南北方向、东南方向三个方向拟合变异函数的套合计算得变程a为2000m，基于变程的方式选取插值数据，平均数据个数为100，这将会引起较大的克里金方程组计算量，甚至方程组不可解。预设地形粗糙度阈值M0为20，此时地形起伏特征明显，基于步长的数据选取方式，水深插值数据95％在步长阈值2(含)以内，此时插值水深数据为12个，运算过程计算量明显降低。

(2)克里金方法的权值修正和克里金插值结果的误差修正

对克里金方程组解得的权值进行修正，从而消除了负权值的水深数据值引起插值结果水深值的偏差；对选取的水深插值数据进行交叉验证，可得水深插值数据的计算值和真值的偏差，用这些水深插值数据的偏差来估计待插值位置水深值的计算偏差，误差修正后的水深插值结果在地形较为平坦区精度相近，在地形变化剧烈区，精度提高明显，从而说明了基于克里金方法的数字高程模型的建立方法的有效性，这对于提取地形匹配辅助导航中的地形特征明显适宜匹配的区域有很好的实用价值。

Claims (1)

1.一种基于克里金方法的水下地形数字高程建立方法，其特征为：

步骤1采用公知的多波束水深测量方法测量水深值，获得水深数据

所述的水深数据表示经度为i、纬度为j的位置A_ij处的水深值，步骤2建立不规则三角网

由初始水深数据在地球坐标系OENU中的NOE平面上的投影点相互连接建立不规则三角网，地球坐标系OE轴正向为地理东向，ON轴正向为地理北向，OU轴正向为天向：

步骤2.1确定投影点的凸壳

选取投影点中纬度值最小的投影点为起点，当存在多个投影点纬度值相等且最小时，取其中经度最小的投影点为起始顶点，记为X₀，由起始顶点X₀向其它投影点且A_ij≠0连线，取其中与X₀连线和OE轴夹角最小的投影点为第一顶点X₁，当存在多个这样的投影点时，取和起始顶点X₀距离值最大的投影点为第一顶点X₁，然后由第一顶点X₁向投影点且A_ij≠0,1连线，取其中与X₁连线和直线X₀X₁夹角最小的投影点为第二顶点X₂,以此类推，选取与X_k的连线和直线X_k-1X_k夹角最小的投影点为第k+1顶点X_k+1，k＝2,3,…n-2，直到所找到的第_n顶点X_n和起始顶点X₀重合为止，_n为所选取得到的顶点个数，由顶点X₀、X₁、X₂…X_n-1依次相连形成凸多边形，并以所述凸多边形为凸壳，顶点X₀、X₁、X₂…X_n-1为凸壳点，

步骤2.2由凸壳建立不规则三角网：

首先，由起始顶点X₀向凸壳的其它顶点依次连线，形成凸壳三角网，然后，对位于凸壳内的投影点作如下连接处理：

任选一个凸壳内的投影点，将所选投影点分别与包围所选投影点的三角形的顶点连线，遍历所有凸壳内的所有投影点，得到由三角形拼接而成的不规则三角网，

步骤2.3根据投影点与水深数据以及由三角形拼接而成的不规则三角网，构建由三角地形单元组成的不规则三角网地形，使得不规则三角网地形与不规则三角网形成对应关系，三角地形单元与不规则三角网中的三角形形成对应关系，步骤3不规则三角网地形空间优化

步骤3.1以任一共边的两个三角地形单元的所有顶点对应的水深数据点作为空间四面体的四个顶点构建空间四面体，所述空间四面体包括共边的两个三角地形单元及新增的两个三角地形单元，

步骤3.2如果空间四面体中的共边的两个三角地形单元的内角标准差之和小于新增的两个三角地形单元的内角标准差之和，则舍弃共边的两个三角地形单元并保留新增的两个三角地形单元；否则，舍弃新增的两个三角地形单元并保留共边的两个三角地形单元，

步骤3.3使用重复步骤3.1和步骤3.2的方式遍历所有共边的两个三角地形单元，得到空间优化的不规则三角网地形，

步骤4采用区域化变量分析方法，对初始水深数据进行结构变异分析，得出变异函数值表达式，所述的区域化变量分析方法为：

步骤4.1分别计算东西向、南北向、东北-西南向各距离的实际变异函数值，

在投影点中，以距离最小的一对投影点的距离值为一个距离单位，距离最大的一对投影点的距离值为h_max，分别筛选出投影点对之间的距离为h′_l＝l×h±△h且过所述投影点对的直线与OE轴夹角θ'＝θ±△θ的所有投影点对，计算筛选出的投影点对所在位置的水深数据对在计算距离为lh的变异函数值，h′_l为距离实际值，h为距离计算值且h＝0.5距离单位～1距离单位，l为距离因子,表示对向下取整运算，△h为距离容限值且△h＝0～0.05h，θ'为角度实际值，θ为角度计算值且东西向、南北向、东北-西南向分别取值为θ＝0°、θ＝90°、θ＝45°，△θ为角度容限值且△θ＝0°～2.5°，实际变异函数值计算如下：

${\mathrm{\&gamma;}}_l=\frac{1}{2N_l}\underset{m=1}{\overset{N_l}{\&Sigma;}}{\&lsqb;Z_{A_{ij}}-Z_{A_{ij}^l}^m\&rsqb;}^2---\left(1\right)$

式中，γ_l为距离值为lh的水深数据对实际变异函数值，N_l为距离值为lh的水深数据对的个数，m为距离值为lh的水深数据对的序号，和分别为位置A_ij处及与其距离为lh的位置处的第m个水深数据对的一个水深值；

步骤4.2拟合实际变异函数值，得到三个方向的变异函数

所述拟合的方法为：选取球型理论变异函数，

${\mathrm{\&gamma;}}_{h^{\&prime;\&prime;}}=\left\{\begin{array}{cc}{C}_0,& h^{\&prime;\&prime;}=0\\ C_0+C(\frac{3}{2}\&CenterDot;\frac{h^{\&prime;\&prime;}}{a}-\frac{1}{2}\&CenterDot;\frac{h^{\&prime;\&prime;3}}{a^3}),& 0<h^{\&prime;\&prime;}\&le;a\\ C_0+C,& h^{\&prime;\&prime;}>a\end{array}\right.---\left(2\right)$

对于东西向、南北向、东北-西南向，计算处对应的理论变异函数值γ_h″，并根据理论变异函数值γ_h″和实际变异函数值γ_l，由最小二乘拟合法，得到三个方向变程值a、块金常数C₀、拱高C，带入球型理论变异函数，得到三个方向的计算变异函数；

步骤4.3三个方向的计算变异函数套合，求取变异函数值计算式，

所述的套合过程如下：

步骤4.3.1坐标系旋转变换

比较东西向、南北向、东北-西南向三个方向的变程值，将原坐标系OE轴和ON轴顺时针旋转β角形成O＇EN系，0°≤β<180°，使OEN坐标系中变程值a最大的方向变为O＇EN系横轴正向，则由OEN系到O＇EN系的旋转矩阵Q的表达式为：

$Q=\left(\begin{array}{cc}cos\mathrm{\&beta;}& sin\mathrm{\&beta;}\\ -sin\mathrm{\&beta;}& \cos \mathrm{\&beta;}\end{array}\right)---\left(3\right)$

步骤4.3.2坐标系压缩变换

将坐标系O＇EN压缩变换成坐标系O＂EN，由O＇EN系到O＂EN系的压缩矩阵P的表达式为：

$P=\left(\begin{array}{cc}1& 0\\ 0& \frac{R_1}{R_2}\end{array}\right)---\left(4\right)$

式中，P为压缩矩阵，R₁,R₂分别为东西向、南北向、东北-西南向三个方向中的最大变程值、最小变程值，

步骤4.3.3求变异函数值计算式

套合后的变异函数值计算式为：

$h^{\&prime;\&prime;\&prime;}=\sqrt{F^T{Q}^T{P}^{T}PQF},F=\left(\begin{array}{c}{i}_0-i_1\\ j_0-j_1\end{array}\right)---\left(5\right)$

${\mathrm{\&gamma;}}_{h^{\&prime;\&prime;\&prime;}}=\left\{\begin{array}{cc}{C_0}^{\&prime;},& h^{\&prime;\&prime;\&prime;}=0\\ {C_0}^{\&prime;}+C^{\&prime;}(\frac{3}{2}\&CenterDot;\frac{h^{\&prime;\&prime;\&prime;}}{R_1}-\frac{1}{2}\&CenterDot;\frac{h^{\&prime;\&prime;\&prime;3}}{{R_1}^3}),& 0<h^{\&prime;\&prime;\&prime;}\&le;R_1\\ {C_0}^{\&prime;}+C^{\&prime;},& h^{\&prime;\&prime;\&prime;}>R_1\end{array}\right.---\left(6\right)$

式中，i₀,j₀,i₁,j₁分别表示在OEN系中的经纬度坐标值，F为由到的向量阵，F^T,P^T,Q^T分别为F,P,Q的转置，h″′为在O＂EN系中与之间的距离，γ_h″′为该距离值h″′的套合后的变异函数值，C₀',C'分别为东西向、南北向、东北-西南向三个变异函数的块金常数、拱高的算术平均值，

步骤5将优化后的不规则三角网地形转换为规则格网地形，转换过程如下：

步骤5.1设定规则格网地形网格间距：

在初始水深数据在OENU系的投影点中，取距离最小的一对投影点的距离值的1/20～1/10为规则格网地形的网格间距，

步骤5.2采用普通克里金插值方法和水深插值误差修正方法估计水深值：

步骤5.2.1采用基于步长的方法，选取普通克里金插值方法的插值数据，所述的基于步长的方法如下：

对于步长的定义为：在不规则三角网中，包围待插值位置的三角形为处的0步长三角形，0步长三角形的顶点为处0步长位置点，与处0步长三角形共边相邻的三角形为处1步长三角形，1步长三角形的非共边顶点为处1步长位置点，以此类推，可得待插值位置处g步长三角形和位置点，和步长三角形相对应的为步长三角地形单元，步长位置点所在的水深数据为步长数据，设定初始粗糙度M₀＝10～30，G为步长阈值，表示步长的最大值，初值为1，计算待插值位置处步长阈值为G＝1内的三角地形单元的地形粗糙度M_G：

$M_G=\frac{\underset{g=0}{\overset{G}{\&Sigma;}}\underset{u_g=1}{\overset{U_g}{\&Sigma;}}{S}_{u_g}}{\underset{g=0}{\overset{G}{\&Sigma;}}\underset{u_g=1}{\overset{U_g}{\&Sigma;}}{S}_{u_g}^{\&prime;}}---\left(7\right)$

式中，u_g为g步长三角形的三角形编号，g＝0,1,…G，U_g为g步长三角形的三角形总数，为第u_g个g步长三角地形单元的空间平面面积，为第u_g个g步长三角形的平面面积，g步长三角形的平面面积为对应的g步长三角地形单元的空间平面面积在OENU系的投影面积，当M_G>M₀,令G＝G+1，进行新的步长阈值G的不规则三角网地形粗糙度M_G的迭代计算，直至迭代结果M_G≤M₀，并选取当前步长阈值G内所有的步长三角地形单元的顶点作为待插值位置处普通克里金插值方法的水深插值数据；

步骤5.2.2普通克里金权系数的计算，

由步骤5.2.1得待插值位置处步长阈值G内的总数为w个水深插值数据根据式(5)，式(6)，计算得待插值位置与水深插值数据在OEN系的投影位置在O＂EN系中的距离值及变异函数值一对水深插值数据在OEN系的投影位置在O＂EN系中的距离值及变异函数值e'＝1,2,…w且e'≠e，带入普通克里金方程组，

$\left[\begin{array}{ccccc}{\mathrm{\&gamma;}}_{h_1^1}& {\mathrm{\&gamma;}}_{h_1^2}& ...& {\mathrm{\&gamma;}}_{h_1^w}& 1\\ {\mathrm{\&gamma;}}_{h_2^1}& {\mathrm{\&gamma;}}_{h_2^2}& ...& {\mathrm{\&gamma;}}_{h_2^w}& 1\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& {\mathrm{\&gamma;}}_{h_e^{e^{\&prime;}}}& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ {\mathrm{\&gamma;}}_{h_w^1}& {\mathrm{\&gamma;}}_{h_w^2}& ...& {\mathrm{\&gamma;}}_{h_w^w}& 1\\ 1& 1& ...& 1& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&lambda;}}_1\\ {\mathrm{\&lambda;}}_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&lambda;}}_w\\ \mathrm{\&mu;}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&gamma;}}_{h_0^1}\\ {\mathrm{\&gamma;}}_{h_0^2}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&gamma;}}_{h_0^w}\\ 1\end{array}\right]---\left(8\right)$

由普通克里金方程组得克里金权系数λ_e,e＝1,2,…w；

步骤5.2.3对普通克里金权系数进行修正，

所述修正的方法为：对权值为负的λ_e置0，其余权值修正为：

${\mathrm{\&lambda;}}_e^{\&prime;}=\frac{{\mathrm{\&lambda;}}_e}{\underset{k=1}{\overset{w}{\&Sigma;}}{\mathrm{\&lambda;}}_k}---\left(9\right)$

式中，λ_e为修正前的权值，λ_k为对负权值进行置0处理后的权值，λ′_e为修正后的权值，

步骤5.2.4普通克里金插值方法的水深插值结果及误差修正，

待插值位置处的水深插值及插值方差的计算：

${\stackrel{\&OverBar;}{Z}}_{A_{ij}^0}=\underset{e=1}{\overset{w}{\&Sigma;}}{\mathrm{\&lambda;}}_e^{\&prime;}{Z}_{A_{ij}^e}---\left(10\right)$

${\mathrm{\&sigma;}}_{A_{ij}^0}=\underset{e=1}{\overset{w}{\&Sigma;}}{\mathrm{\&lambda;}}_e^{\&prime;}{\mathrm{\&gamma;}}_{h_0^e}+\mathrm{\&mu;}---\left(11\right)$

对待插值位置处每一个水深插值数据进行水深值交叉验证，取水深插值数据由式(5)、式(6)、式(7)、式(8)、式(9)、式(10)、式(11)计算得位置处的水深计算值及误差

修正后位置处的普通克里金插值水深值为：

${\stackrel{\&OverBar;}{Z}}_{A_{ij}^0}^{\&prime;}={\stackrel{\&OverBar;}{Z}}_{A_{ij}^0}+\frac{{\mathrm{\&sigma;}}_{A_{ij}^0}}{\underset{e=1}{\overset{w}{\&Sigma;}}{\mathrm{\&sigma;}}_{A_{ij}^e}}\&CenterDot;\underset{e=1}{\overset{w}{\&Sigma;}}(Z_{A_{ij}^e}-{\stackrel{\&OverBar;}{Z}}_{A_{ij}^e})---\left(12\right)$

步骤6输出规则格网数字地形。