CN106484997A - 一种基于克里金插值的水岸带淤泥厚度计算及出图方法 - Google Patents

Abstract

本发明公开了一种基于克里金插值的水岸带淤泥厚度计算及出图方法：首先使用船载RTK和声呐进行点数据采集；然后，分别对高频和低频点数据进行普通克里金格网插值，形成两个频段的规则格网数据；最后，结合杆高等固定参数设置，计算淤泥顶和淤泥底数据，得到淤泥厚度，并使用ArcGIS进行栅格计算，得到空间淤泥厚度图。本发明的基于克里金插值的水岸带淤泥厚度计算及出图方法，将无人船应用于河道的淤泥厚度测绘中，能够提高测绘的工作效率和工作准确度。通过对南京水西门河段进行实地测试和出图验证，本方法具有较高的精度和效率，能够较好的完成典型水岸带河底淤泥的精确测定，能够为水岸带土壤研究提供技术支持。

Description

一种基于克里金插值的水岸带淤泥厚度计算及出图方法

技术领域

本发明属于测量淤泥厚度技术领域，具体涉及一种基于克里金插值的水岸带淤泥厚度计算及出图方法。

背景技术

湖滨水岸带是介于陆生生态系统与水生生态系统之间的生态交错类型之一，生产力高，生态边缘效应显著。湖滨水岸带河底土壤长期处于厌氧与好氧交替的环境中，强烈积累有机物质并且分解缓慢，故有机质丰富。对水岸带河底淤泥的厚度测定对研究水岸带土壤特性及活性有机碳采样提取有着重要的参考意义。原始的淤泥厚度测量方法，一般情况下操作人员必须下水，直接面对河水及杂草、河底淤泥，工作环境比较恶劣，承载的测量船只也有倾覆的危险。

近年来，测绘无人船技术得到了较大发展，无人船不仅能够对搭载水质监测仪，对水质进行实时监测，也能够对河底进行深度探测。而利用测绘无人船对于淤泥厚度的研究，由于对设备要求较高，目前工作中还缺乏实际的操作经验。

发明内容

发明目的：针对现有技术中存在的不足，本发明的目的是提供一种基于克里金插值的水岸带淤泥厚度计算及出图方法，提高测绘的工作效率和工作准确度。

技术方案：为了实现上述发明目的，本发明采用的技术方案为：

一种基于克里金插值的水岸带淤泥厚度计算及出图方法：首先使用船载RTK和声呐进行点数据采集；然后，分别对高频和低频点数据进行普通克里金格网插值，形成两个频段的规则格网数据；最后，结合杆高等固定参数设置，计算淤泥顶和淤泥底数据，得到淤泥厚度，并使用ArcGIS进行栅格计算，得到空间淤泥厚度图。

所述的基于克里金插值的水岸带淤泥厚度计算及出图方法，包括以下步骤：

1)通过岸端服务器进行测绘无人船进行路线及目标点设置；确定河道的起始点和终止点，然后设置无人船测量的时间间隔为1次/秒，采样路线采用Z字形采样路线；将预先设定好的路线输入控制系统，让无人船能够进行自动巡航和测量数据，所测量得到的数据是以CAD软件的DXF格式保存；

2)无人船将所测量的CAD的DXF点状数据，通过无线传输的方式，传输到岸端服务器控制电脑中，利用ArcGIS软件的格式转换工具，转为shp格式；

3)设置克里金参数并插值，生成格网数据；

4)利用克里金插值计算得到的高频和低频信号的栅格数据，结合杆高、水面高、无人船吃水深度固定参数设置，计算淤泥顶部和底部数据；

5)在获得淤泥底高和顶高的基础上，将顶高减去底高，就得到了淤泥的厚度图；

6)最后，将实际人工测量数据与所得到的淤泥厚度进行验证，并使用ArcGIS软件进行出图。

步骤3)中，使用克里金插值将测量点构建成测量面；将所要测试的内河河道段区域设定为B，将需要构面的河道淤泥层测量值设定为Z(x)，因此有{Z(x)∈B}，Z(x)是一个二阶平稳的随机函数，它在空间中的取值设定为Z(x₁),Z(x₂),...,Z(x_n)，其中x表示淤泥测量点的空间位置；按照普通克里金算法的原理，非测量点x₀的高频和低频值的Z(x₀)估计值是多个已知测量点的加权求和，如下式所示：

其中，Z(x_i)(i＝0,1,2,...,n)为测量点x_i的值，x₀为非测量点，也就是需要插值得到的淤泥层值，其他为已知点；λ_i(i＝1,2,...,n)为权值；权值并非仅由距离决定，是需要在最小方差及无偏特征的假设条件下，由变异函数计算数值决定；

普通克里金方程组：

按照普通克里金估计方差最小，可以变形成为如下所示：

其中，C表示协方差函数，E{}为数学期望值。

步骤4)中，使用的克里金插值算法进行淤泥测量点计算，步骤具体如下：

1)对于观测数据，两两计算距离与半方差；

2)寻找一个拟合曲线拟合距离与半方差的关系，从而能根据任意距离计算出相应的半方差；

3)计算出所有已知点之间的半方差r_ij；

4)对于未知点z_o，计算它到所有已知点z_i的半方差r_io；

5)求解步骤4)的方程组，得到最优系数λi；

6)使用最优系数对已知点的属性值进行加权求和，得到未知点z_o的估计值。

步骤5)中，具体公式如下所示：

H_淤泥顶高＝(H_水面高程-H_高频-(H_杆高+H_吃水高度))

H_淤泥底高＝(H_水面高程-H_低频-(H_杆高+H_吃水高度))

H_淤泥厚度＝H_淤泥顶高-H_淤泥底高。

有益效果：与现有技术相比，本发明的基于克里金插值的水岸带淤泥厚度计算及出图方法，将无人船应用于河道的淤泥厚度测绘中，能够提高测绘的工作效率和工作准确度。通过对南京水西门河段进行实地测试和出图验证，本方法具有较高的精度和效率，能够较好的完成典型水岸带河底淤泥的精确测定，能够为水岸带土壤研究提供技术支持。

附图说明

图1是基于克里金插值的水岸带淤泥厚度计算及出图方法的流程图；

图2是克里金插值结果图；其中，a为高频面结果图，b为低频面结果图；亮度越暗，表示深度越深；

图3是河道淤泥厚度图。

具体实施方式

下面结合具体实施例对本发明做进一步的说明。

以下实施例中所使用的测绘无人船，长约1.5米，宽0.7米，船身安装了RTK装置及双频声呐装置。

RTK(Real Time Kinematic，载波相位差分)定位技术是基于载波相位观测值的实时动态定位技术，能够实时地提供测站点在指定坐标系中的三维定位结果，并达到厘米级精度。精度高，作业方便。RTK作业不受通视条件限制，无需做控制，基准站设置好，进行点检核后，即可开测，如用虚拟基站则更简便。由于精度达到厘米级，所以可以用于河道淤泥的测量点的精确定位。

双频声呐，指的是能够在同一时段同一地点发射2个频段的声呐信号，利用声波的频率和强度存在的差异，测定所要求的被测物体的形态和距离等参数特征。由于声波信号的频率不同，其在水下的穿透力和强度也存在差异。高频声呐信号设定为180MHz，穿透力弱，信号强，能够穿透内陆河道中的浮游植物和悬浮垃圾等，但是无法穿透淤泥顶部，可以用来测定声呐距离淤泥顶部的长度；低频声呐信号设定为50MHz，穿透力较强，信号偏弱，它不仅能够轻易穿透浮游植物和悬浮垃圾，也能够穿透淤泥层，直接射到淤泥底部的河床硬质层。由于河床硬质层的密度大，低频声呐信号无法穿透，因此可以用来测定声呐距离河床硬质层的长度，也就是距离淤泥底部的长度。

实施例1

一种基于克里金插值的水岸带淤泥厚度计算及出图方法，流程图如图1的所示，具体步骤如下：

1)测量淤泥样本点数据，生成CAD数据。

通过岸端服务器进行测绘无人船进行路线及目标点设置。确定河道的起始点和终止点，然后设置无人船测量的时间间隔为1次/秒，采样路线采用Z字形采样路线。将预先设定好的路线输入控制系统，让无人船能够进行自动巡航和测量数据，所测量得到的数据是以CAD软件的DXF格式保存。

2)将CAD数据转换成shp点状格式数据。

shp格式是shapefile格式的简写，该格式是由ESRI公司开发，一个shp文件包括一个主文件，一个索引文件，和一个dBASE表。之所以要将原始格式DXF转换成shp格式，是因为ArcGIS软件具有强大的出图功能，最后出图将会使用ArcGIS软件进行操作，而ArcGIS软件无法直接操作DXF文件，但是能够直接操作shp格式数据。无人船将所测量的CAD的DXF点状数据，通过无线传输的方式，传输到岸端服务器控制电脑中，利用ArcGIS软件的格式转换工具，转为shp格式。

3)设置克里金参数并插值，生成格网数据。

由于克里金插值能够将点状数据生成栅格格网数据，并且能够对每个数据区域的精度进行求解，所以采用克里金插值将点生成面。所使用的软件工具是ArcGIS软件，计算得到生成高频和低频规则格网的栅格数据，用于后期的求取层差，从而获得淤泥厚度。

4)基于格网数据，计算淤泥顶高和淤泥底高的格网数据。

利用克里金插值计算得到的高频和低频信号的栅格数据，结合杆高、水面高、无人船吃水深度等固定参数设置，计算淤泥顶部和底部数据。

5)基于淤泥底高和顶高的格网数据，进行差值计算，获得淤泥厚度图。

在获得淤泥底高和顶高的基础上，将顶高减去底高，就得到了淤泥的厚度图。采用的具体操作是实用栅格相减。

6)进行结果验证，利用ArcGIS软件进行出图。

最后，将实际人工测量数据与所得到的淤泥厚度进行验证，并使用ArcGIS软件进行出图。

本专利方法使用克里金(Kriging)插值将测量点构建成测量面。克里金插值首先考虑的是空间属性在二维空间中的分布，然后用此范围内的采样点来估计待插点的属性值。将所要测试的内河河道段区域设定为B，将需要构面的河道淤泥层测量值设定为Z(x)，因此有{Z(x)∈B}，Z(x)是一个二阶平稳的随机函数，它在空间中的取值设定为Z(x₁),Z(x₂),...,Z(x_n)，其中x表示淤泥测量点的空间位置。按照普通克里金算法的原理，非测量点x₀的高频和低频值的Z(x₀)估计值是多个已知测量点的加权求和，如下式所示：

其中，Z(x_i)(i＝0,1,2,...,n)为测量点x_i的值，x₀为非测量点，也就是我们需要插值得到的淤泥层值，其他为已知点。λ_i(i＝1,2,...,n)为权值。权值并非仅由距离决定，是需要在最小方差及无偏特征的假设条件下，由变异函数计算数值决定。

详细推导请参见文献[靳国栋,刘衍聪,牛文杰.距离加权反比插值法和克里金插值法的比较[J].长春工业大学学报,2003,24(3):53-57]。以下列出普通克里金方程组：

按照普通克里金估计方差最小，可以变形成为如下所示：

其中，C表示协方差函数，E{}为数学期望值。

本实施例使用的普通克里金插值算法进行淤泥测量点计算的基本步骤，具体如下：

1)对于观测数据，两两计算距离与半方差；

2)寻找一个拟合曲线拟合距离与半方差的关系，从而能根据任意距离计算出相应的半方差；

3)计算出所有已知点之间的半方差r_ij；

4)对于未知点z_o，计算它到所有已知点z_i的半方差r_io；

5)求解步骤4)的方程组，得到最优系数λi；

6)使用最优系数对已知点的属性值进行加权求和，得到未知点z_o的估计值。

本专利将游离于河床表面，一定密度内的底质层称之为河道淤泥层。通常其组成为粘性细颗粒泥沙。本实施例使用的双频声呐采用的工作原理是：声呐装置同时使用2个工作频率，180kHz和50kHz。其中，低频由于其波长长，穿透力强，能够到达河床底；高频的穿透能力较弱，只能到达淤泥顶部。然后，利用两者的高程差，就能够计算出淤泥的厚度。具体公式如下所示：

H_淤泥顶高＝(H_水面高程-H_高频-(H_杆高+H_吃水高度))

H_淤泥底高＝(H_水面高程-H_低频-(H_杆高+H_吃水高度))

H_淤泥厚度＝H_淤泥顶高-H_淤泥底高

实施例2

以江苏省南京市水西门河段为例，采用实施例1的方法进行实地测量，本实施例使用的计算机硬件主要为：CPU是i5处理器，内存是8GB，硬盘为固态硬盘。使用的软件为：操作系统是windows 7专业版系统，ArcGIS软件是9.3.1版本。

(1)首先是对采集到的点状数据进行格式转化。为了使用ArcGIS软件中的克里金插值算法，需要将原始格式转化为shapefile格式数据。测绘无人船直接测量的数据格式是CAD的DXF格式，使用ArcGIS软件将其打开后，选中“*.dxf Point”后，右键点击，在弹出窗口中选择Data，然后是Export Data，就可以将其转化为shapefile格式。

(2)其次是构建克里金规则格网栅格数据。由于高频点和低频点都是点格式，无法直接利用，也无法直接生成面数据，这就需要使用克里金插值算法构建规则格网栅格面。使用ArcGIS软件中的Geostatistical Analyst(地统计分析工具)进行计算，采用的是经典的普通克里金插值算法。插值后的高频面和低频面见图2所示，其中，a为高频面结果图，b为低频面结果图；亮度越暗，表示深度越深。

(3)使用高频面和低频面的插值，计算淤泥厚度图。如果需要计算淤泥顶高和淤泥底高，则需要知道参数杆高和吃水深度。由于只需要获得淤泥厚度，因此无需知道杆高和吃水深度，直接用高频面减去低频面，就能得到淤泥厚度。使用ArcGIS软件中的RasterCalculator进行运算，得到的结果如图3所示。使用10级分级图显示，亮度越高，则淤泥厚度约薄，亮度越暗，则淤泥厚度越厚。从结果图可知，河道的两侧的淤泥厚度更厚，而河道中央的淤泥厚度则相对较薄，这是由于南京水西门的外秦淮河的水流较为湍急，中间的淤泥都被冲刷到河道的两侧，这就为今后该河段的清淤工作指明了方向，清理淤泥需要从河道两侧重点处理。

表1所示的是低频面、高频面以及淤泥厚度面的统计值。从中可以发现，低频面可以近似指代淤泥底部，均值为3.58米；高频面可以近似指代淤泥顶部，均值为7.76米。淤泥顶部的标准差为0.01，说明其淤泥分布较为均匀。

表1实验区淤泥面统计值

可见，使用测绘无人船进行河底淤泥厚度的测量，不仅能够减少人力物力的投资，降低测量人员的测绘风险，也能够提高测绘的效率和精度，是未来的发展趋势。从测量的结果可知，本实验所测的水西门地区外秦淮河河道的两侧的淤泥厚度更厚，而河道中央的淤泥厚度则相对较薄。

Claims (5)

1.一种基于克里金插值的水岸带淤泥厚度计算及出图方法，其特征在于：首先使用船载RTK和声呐进行点数据采集；然后，分别对高频和低频点数据进行普通克里金格网插值，形成两个频段的规则格网数据；最后，结合杆高等固定参数设置，计算淤泥顶和淤泥底数据，得到淤泥厚度，并使用ArcGIS进行栅格计算，得到空间淤泥厚度图。

2.根据权利要求1所述的基于克里金插值的水岸带淤泥厚度计算及出图方法，其特征在于：包括以下步骤：

1)通过岸端服务器进行测绘无人船进行路线及目标点设置；确定河道的起始点和终止点，然后设置无人船测量的时间间隔为1次/秒，采样路线采用Z字形采样路线；将预先设定好的路线输入控制系统，让无人船能够进行自动巡航和测量数据，所测量得到的数据是以CAD软件的DXF格式保存；

2)无人船将所测量的CAD的DXF点状数据，通过无线传输的方式，传输到岸端服务器控制电脑中，利用ArcGIS软件的格式转换工具，转为shp格式；

3)设置克里金参数并插值，生成格网数据；

4)利用克里金插值计算得到的高频和低频信号的栅格数据，结合杆高、水面高、无人船吃水深度固定参数设置，计算淤泥顶部和底部数据；

5)在获得淤泥底高和顶高的基础上，将顶高减去底高，就得到了淤泥的厚度图；

6)最后，将实际人工测量数据与所得到的淤泥厚度进行验证，并使用ArcGIS软件进行出图。

3.根据权利要求2所述的基于克里金插值的水岸带淤泥厚度计算及出图方法，其特征在于：步骤3)中，使用克里金插值将测量点构建成测量面；将所要测试的内河河道段区域设定为B，将需要构面的河道淤泥层测量值设定为Z(x)，因此有{Z(x)∈B}，Z(x)是一个二阶平稳的随机函数，它在空间中的取值设定为Z(x₁),Z(x₂),...,Z(x_n)，其中x表示淤泥测量点的空间位置；按照普通克里金算法的原理，非测量点x₀的高频和低频值的Z(x₀)估计值是多个已知测量点的加权求和，如下式所示：

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_{i}Z\left(x_i\right)$

其中，Z(x_i)(i＝0,1,2,...,n)为测量点x_i的值，x₀为非测量点，也就是需要插值得到的淤泥层值，其他为已知点；λ_i(i＝1,2,...,n)为权值；权值并非仅由距离决定，是需要在最小方差及无偏特征的假设条件下，由变异函数计算数值决定；

普通克里金方程组：

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_{i}C(x_i,x_j)-u=C(x_0,x_j),j=1,2,...,n\\ \underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i=1\end{array}\right.$

按照普通克里金估计方差最小，可以变形成为如下所示：

${\mathrm{\σ}}_{OK}^2=E\left\{{\[Z\left(x_0\right)-Z^{*}\left(x_0\right)\]}^2\right\}=C(x_0,x_j)-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i}C(x_0,x_i)+u$

其中，C表示协方差函数，E{}为数学期望值。

4.根据权利要求2所述的基于克里金插值的水岸带淤泥厚度计算及出图方法，其特征在于：步骤4)中，使用的克里金插值算法进行淤泥测量点计算，步骤具体如下：

1)对于观测数据，两两计算距离与半方差；

2)寻找一个拟合曲线拟合距离与半方差的关系，从而能根据任意距离计算出相应的半方差；

3)计算出所有已知点之间的半方差r_ij；

4)对于未知点z_o，计算它到所有已知点z_i的半方差r_io；

5)求解步骤4)的方程组，得到最优系数λi；

6)使用最优系数对已知点的属性值进行加权求和，得到未知点z_o的估计值。

5.根据权利要求2所述的基于克里金插值的水岸带淤泥厚度计算及出图方法，其特征在于：步骤5)中，具体公式如下所示：

H_淤泥顶高＝(H_水面高程-H_高频-(H_杆高+H_吃水高度))

H_淤泥底高＝(H_水面高程-H_低频-(H_杆高+H_吃水高度))

H_淤泥厚度＝H_淤泥顶高-H_淤泥底高。