[发明专利]一种基于克里金方法的水下地形数字高程建立方法

摘要

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

主权项

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遍历所有共边的两个三角地形单元，得到空间优化的不规则三角网地形，步骤4采用区域化变量分析方法，对初始水深数据进行结构变异分析，得出变异函数值表达式，所述的区域化变量分析方法为：步骤4.1分别计算东西向、南北向、东北-西南向各距离的实际变异函数值，在投影点中，以距离最小的一对投影点的距离值为一个距离单位，距离最大的一对投影点的距离值为h_max，分别筛选出投影点对之间的距离为h′_h＝l×h±Δh且过所述投影点对的直线与OE轴夹角θ′＝θ±Δθ的所有投影点对，计算筛选出的投影点对所在位置的水深数据对在计算距离为lh的变异函数值，h′_l为距离实际值，h为距离计算值且h＝0.5距离单位～1距离单位，l为距离因子，表示对向下取整运算，△h为距离容限值且Δh＝0～0.05h，θ′为角度实际值，θ为角度计算值且东西向、南北向、东北-西南向分别取值为θ＝0°、θ＝90°、θ＝45°，Δθ为角度容限值且Δθ＝0°～2.5°，实际变异函数值计算如下：γl=12NlΣm=1Nl[ZAij-ZAijlm]2---(1)]]>式中，γ_l为距离值为lh的水深数据对实际变异函数值，N_l为距离值为lh的水深数据对的个数，m为距离值为lh的水深数据对的序号，和分别为位置A_ij处及与其距离为lh的位置处的第m个水深数据对的一个水深值；步骤4.2拟合实际变异函数值，得到三个方向的变异函数所述拟合的方法为：选取球型理论变异函数，γh′′=C0,h′′=0C0+C(32·h′′a-12·h′′3a3),0<h′′≤aC0+C,h′′>a---(2)]]>对于东西向、南北向、东北-西南向，计算处对应的理论变异函数值γ_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=cosβsinβ-sinβcosβ---(3)]]>步骤4.3.2坐标系压缩变换将坐标系O′EN压缩变换成坐标系O″EN，由O′EN系到O″EN系的压缩矩阵P的表达式为：P=100R1R2---(4)]]>式中，P为压缩矩阵，R₁，R₂分别为东西向、南北向、东北-西南向三个方向中的最大变程值、最小变程值，步骤4.3.3求变异函数值计算式套合后的变异函数值计算式为：h′′′=FTQTPTPQF,F=i0-i1j0-j1---(5)]]>γh′′′=C0′,h′′′=0C0′+C′(32·h′′′R1-12·h′′′3R13),0<h′′′≤R1C0′+C′,h′′′>R1---(6)]]>式中，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：MG=Σg=0GΣug=1UgSugΣg=0GΣug=1UgSug′---(7)]]>式中，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的不规则三角网地形粗糙度MG的迭代计算，直至迭代结果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，带入普通克里金方程组，γh11γh12...γh1w1γh21γh22...γh2w1.........γhee′.........γhw1γhw2...γhww111...10λ1λ2...λwμ=γh01γh02...γh0w1---(8)]]>由普通克里金方程组得克里金权系数λ_e，e＝1，2，…w；步骤5.2.3对普通克里金权系数进行修正，所述修正的方法为：对权值为负的λ_e置0，其余权值修正为：λe′=λeΣk=1wλk---(9)]]>式中，λ_e为修正前的权值，λ_k为对负权值进行置0处理后的权值，λ′_e为修正后的权值，步骤5.2.4普通克里金插值方法的水深插值结果及误差修正，待插值位置处的水深插值及插值方差的计算：Z‾Aij0=Σe=1wλe′ZAije---(10)]]>σAij0=Σe=1wλe′γh0e+μ---(11)]]>对待插值位置处每一个水深插值数据e＝1，2，…w进行水深值交叉验证，取水深插值数据由式(5)、式(6)、式(7)、式(8)、式(9)、式(10)、式(11)计算得位置处的水深计算值及误差修正后位置处的普通克里金插值水深值为：Z‾Aij0′=Z‾Aij0+σAij0Σe=1wσAije·Σe=1w(ZAije-Z‾Aije)---(12)]]>步骤6输出规则格网数字地形。