多波束测深数据处理的抗差最小二乘配置迭代解法

2017-06-07王乐洋陈汉清

测绘学报 2017年5期

王乐洋，陈汉清

1. 东华理工大学测绘工程学院，江西南昌 330013； 2. 流域生态与地理环境监测国家测绘地理信息局重点实验室，江西南昌 330013； 3. 江西省数字国土重点实验室，江西南昌330013

王乐洋1,2,3，陈汉清1,2

针对利用最小二乘配置处理多波束测深数据，存在二次曲面数学模型通常无法精确表征海底地形的整体变化趋势以及观测数据存在粗差或异常点时，常规方法给出的协方差函数不能精确表征其统计特性的问题，本文提出了一种抗差最小二乘配置迭代解法。该方法首先进行协方差函数和观测值方差阵初始化，以多面函数拟合趋势项，然后应用等价权抗差估计并通过迭代计算，最终给出稳健的协方差函数参数解及最小二乘配置解。利用本文提出的方法及传统的方法处理实测的多波束测深数据，试验结果表明，相比于传统的方法，本文提出的方法能够较好地表征海底地形的整体变化趋势，一定程度上克服了多波束测深数据中粗差或异常点的影响。相比于传统的抗差方法，本文方法更为有效地识别出测深数据中异常点，推估效果较好，具有稳健性。

最小二乘配置；海底地形生成；协方差函数；抗差；趋势项

近年来，在海底地形测量技术中多波束测深系统成为海洋测量的主要设备，给研究区域的海底地形构建带来了丰富的数据资料[1]。高精度的海底地形地貌信息是保障海上船只航行及舰艇作战训练安全的必备资料，也是海洋工程、资源开发和海洋科学研究等多个应用领域的重要基础信息[2]。为了提高海底地形构建的精度，必须采用合适的算法对多波束测深数据进行处理，众多学者作了大量的研究工作[3]。文献[3]提出了抗差Kriging拟合方法，该方法在精确估计的同时能够检测测深数据的异常点，但是该方法需要测深数据满足平稳假设或内蕴假设，只能对平坦的海底数据进行处理。

最小二乘配置最初是由Trarup提出，是用来研究地球形状和重力场的一种数学方法[4]，随后Moritz和Wolf对最小二乘配置作了进一步研究[5-6]。如今，最小二乘配置已被广泛应用于地壳形变分析[7-9]、高程异常[10]、重力异常[11]、坐标转换[12]等领域。文献[13]将最小二乘配置应用于海道测量非定位点的估计中，利用区域深度场信息解算出非定位点的水深值，该方法方便精度评估，取得了不错的效果。最小二乘配置能够顾及多波束测深数据具有趋势性和随机性的影响，但无法探测出多波束测深数据中的粗差或异常点，针对于此问题，文献[14]提出了一种抗差最小二乘配置方法进行处理，取得了良好的效果。文献[15]对多波束测深数据处理的传统抗差最小二乘配置方法和抗差Kriging拟合方法进行了比较分析，通过试验验证了在多波束测深数据处理中，最小二乘配置在模型上优于Kriging拟合方法。但现有的多波束测深数据处理的配置法均基于二次多项式拟合趋势项，一般情况下，特别是在海底地形起伏较大时，二次多项式难以表征海底地形的整体变化趋势，甚至拟合的地形可能会出现“龙格现象”[16]，从而导致后续的经验协方差函数估计存在偏差；另一方面，由于海况环境、测量仪器以及人为等因素的影响，多波束测深数据不可避免的存在粗差点或异常点的情况，迫切需要一种有效的抗差方法进行处理。

为了解决以上两个问题，提高多波束测深数据处理的精度，本文提出一种抗差最小二乘配置迭代解法，该方法以多面函数代替传统的二次多项式拟合海底地形的趋势项，应用Huber权函数并通过迭代的方式求出经验协方差函数参数估计和最小二乘配置解。最后利用本文方法、传统的抗差最小二乘配置方法及最小二乘配置方法进行实际的多波束测深数据处理，验证了本文方法的有效性与可行性。

1 多波束测深数据处理的最小二乘配置方法

最小二乘配置方法能够同时顾及到海底地形的趋势项及随机项，根据海底地形生成的特点，建立多波束测深观测模型[17-18]

Z=GY+BX+Δ

(1)

(2)

式中，DΔ为测深值方差阵，PΔ为测深值权阵。将式(1)改写成误差方程形式

(3)

根据广义最小二乘原则，构造目标函数

(4)

式中，PX为多波束测深已测点与未测点的权阵，由式(4)求自由极值，并通过整理，得

(5)

将式(3)代入式(5)，通过法方程解算，得到多波束测深数据处理的最小二乘配置解[17]

(6)

式中，P=(DΔ+DS)-1；DS为测深值已测点信号协方差；DS′S为测深值未测点信号与已测点信号的协方差。其中，DS和DS′S由经验协方差函数给出，常选定高斯函数进行拟合[8-9,14]

(7)

(8)

(9)

式中，(xi,yi)为i点测深值的坐标，Y=[abc

def]T为二次多项式模型待估参数。但一般情况下，特别是在海底地形起伏较大的区域，利用二次多项式拟合海底的整体连续变化趋势可能会出现“龙格现象”[16]，不能够精确地表征海底地形的整体连续变化情况。本文针对此问题，提出利用多面函数代替传统的二次多项式进行海底地形整体变化趋势的拟合。

文献[19—20]提出了多面函数模型来逼近不规则表面的思想，即任何的一个数学表面均可由一系列的数学表面的总和，以任意的精度进行逼近。文献[21]提出了基于多面函数的最小二乘配置改进方法，该方法利用多面函数拟合地壳形变的趋势项，并得到了较好的效果。假设海底地形的整体趋势是一个连续不规则变化的表面，基于多面函数的优良特性，利用多面函数对海底地形的整体趋势项进行逼近[20]

(10)

式中，γi为待定系数，R为核函数，(xi,yi)为核函数的选定点坐标。其中，核函数R的表达形式为

(11)

h=Rγ

(12)

以多面函数拟合海底地形的整体连续变化趋势，在最小二乘配置的框架下，把多面函数融入最小二乘配置中，此时基于多面函数的最小二乘配置函数模型为

Z=Rγ+BX+Δ

(13)

式中，Rγ表示海底地形的整体连续变化部分；R为n×m维的核函数矩阵；γ为m×1维的待定系数向量；BX为海底地形不规则变化的随机部分。同理，根据广义最小二乘原理，解得基于多面函数的最小二乘配置解[21]

(14)

2 多波束测深数据处理的抗差最小二乘配置迭代解法

2.1 观测值初始权阵的确定

多波束测深系统是现今海洋测深的主要数据采集设备，它获取的数据覆盖范围广，采样密度高等特点而被广泛应用。但由于海洋测量相比于陆地测量存在较多的不稳定因素，导致多波束测深数据存在粗差或异常点的情况较多，当多波束测深数据含有粗差或异常点时，给出的最小二乘配置解并非是一个最优解。针对多波束测深数据含有粗差或异常点的问题，在文献[21—24]的基础上，本文提出了抗差最小二乘配置方法处理含有粗差或异常点的多波束测深数据。该方法首先进行观测值权阵的初始化，具体权阵初始化步骤如下：

(1) 对多波束测深数据消除整体连续变化趋势，求取初始残差值

(15)

(2) 根据给出的残差初始值，计算初始方差因子，利用该式去计算参数的抗差解具有50%的高崩溃污染率

(16)

(3) 选取Huber权函数计算权因子[14]

(17)

2.2 协方差函数的抗差拟合方法

经验协方差函数拟合的准确性是求解最小二乘配置问题的关键。经验协方差函数的拟合主要为求解协方差函数(式(7))的两个参数解D(0)和k。由上述得到的更新后的观测值权阵，代入式(8)，得到抗差后的D(0)解

(18)

(19)

(20)

(1) 首先对式(20)左边的k2开根号处理，并取其中位数，计算其残差值，即

(21)

(2) 根据给出的残差初始值，计算初始方差因子

(22)

(3) 选取Huber权函数计算权因子

(23)

最后通过迭代计算，给出稳健的协方差函数参数k的值。抗差最小二乘配置迭代解法的计算步骤如下：

(1) 初始化协方差函数，给出协方差函数的初始值D(0)(0)=1，k(0)=0。利用观测值中位数，初始观测值方差阵，以多面函数拟合趋势项，计算得到最小二乘配置的初始解；

(2) 由步骤(1)得到的初始最小二乘配置解，去除趋势项，根据Huber权函数更新观测值方差阵，随后进行信号协方差统计，根据3.2节的方法得到协方差函数的参数解D(0)和k；

(3) 由步骤(2)得到的新的协方差函数以及更新后的观测值方差阵，计算信号协方差Ds，把Ds和更新后的观测值方差阵代入式(14)，计算得到新的最小二乘配置解。

3 算例分析

本文数据来自于中国某海域，水深范围为1100～1400 m，海底地形起伏变化较大。本文选取其中一个条带作为试验数据，由于数据量较大，冗余点较多，首先进行了数据抽稀处理，共选取886个多波束测深数据点，其中随机抽取不含粗差的150个点作为检核点，所有数据如图1所示，星点为检核点，其余为已测点。测深值等高线图如图2所示。

图1 多波束测深数据分布图Fig.1 Multi-beam data distribution map

图2 原始数据测深等高线图Fig.2 Contour map of raw data

为了验证本文算法的有效性，分别利用本文算法、传统抗差最小二乘配置方法、基于多面函数的最小二乘配置方法和基于二次多项式的最小二乘配置方法进行多波束测深数据处理，设计的4种方案如下所示：

方案1：基于二次多项式的最小二乘配置法，该法以二次多项式拟合趋势项，计算检核点的水深值和外符精度；

方案2：基于二次多项式的传统抗差最小二乘配置法，计算检核点的水深值和外符精度；

方案3：基于多面函数的最小二乘配置法，该法以多面函数拟合趋势项，计算检核点的水深值和外符精度；

方案4：基于多面函数的抗差最小二乘配置迭代解法，该法以多面的函数拟合趋势面，结合抗差的方法进行协方差函数拟合，计算检核点的水深值和外符精度。

通过以上4种方案对实际多波束测深数据进行处理，4种方案的外符精度如表1所示，残差统计如图3所示。

表1 4种方案检核点的残差统计情况

Tab.1 Residual statistics of external check points using four schemes m

通过表1可以看出，基于多面函数的最小二乘配置方法的各项指标均比传统的基于二次多项式的最小二乘配置要优。在外符精度上，基于多面函数的最小二乘配置(方案3)和基于二次多项式的最小二乘配置方法(方案1)的外符精度分别为0.935 86 m和1.205 97 m，相比于传统的方法，基于多面函数的最小二乘配置方法在外符精度上提高了0.270 11 m，说明了以多面函数拟合海底地形的趋势项更能反映研究区域海底地形的整体连续变化的趋势，去除趋势项后得到的随机变化部分也更为平稳，因此可以利用该方法得到更为正确的信号值。通过比较两种方法的残差统计图(图3)，可以看出基于多面函数的最小二乘配置方法给出的残差小于传统的基于二次多项式的最小二乘配置方法的残差，进一步验证多面函数拟合趋势项的有效性。另一方面，通过比较分析传统的抗差最小二乘配置方法(方案2)与基于多面函数的最小二乘配置方法(方案3)可以发现，虽然传统的抗差方法进行了抗差处理，但是由于二次多项式不能较好地表征海底地形的整个连续变化趋势，使得后续由扣除趋势项后得到的残差进行拟合得到的协方差函数发生了偏差，导致推估的效果比没有考虑到抗差的基于多面函数的最小二乘配置方法要差。最后，本文抗差最小二乘配置方法的推估效果最好，在外符精度上，比基于二次多项式的最小二乘配置方法提高了0.308 13 m，比传统的抗差最小二乘配置方法提高了0.317 94 m，比基于多面函数的最小二乘配置方法提高了0.038 02 m，说明了本文抗差方法能在一定程度上克服多波束测深数据中粗差或异常点对最小二乘配置结果的影响；同时，相比传统的抗差方法，本文方法推估精度更高。同样的，为了进一步验证本文抗差方法的有效性的，给出四种方案的经验协方差函数的参数解，具体如表2所示。本文抗差最小二乘配置能够在一定程度上抵抗多波束测深数据中粗差或异常点的影响，给出协方差函数更符合实际情况，更为客观的反映点位之间的相关性。

表2 4种方案的协方差函数参数解

Tab.2 Fitting results of covariance function parameters using four schemes

参数方案1方案2方案3方案4D(0)13.9217813.006262.440060.93110k26.95240×10-86.87077×10-85.97079×10-83.35709×10-7

图3 4种方案检核点残差统计图Fig.3 Residual statistics chart of external check points using four schemes

图4 原始多波束测深数据TIN图Fig.4 TIN chart of original multi-beam data

图5 传统抗差方法处理后的多波束测深数据TIN图Fig.5 TIN chart of multi-beam data using traditional robust least squares collocation

图6 经本文抗差方法处理后的TIN图Fig.6 TIN chart of multi-beam data using robust least squares collocation was proposed in this paper

4 结论

为了进一步提高多波束测深数据处理的精度，本文提出了一种抗差最小二乘配置迭代解法，给出了具体的求解公式以及迭代算法，通过算例验证了本文方法的有效性，得出了以下的结论：

(1) 针对现有配置法的趋势项的二次曲面数学模型通常无法精确表征海底地形的整体连续变化趋势的问题，本文提出以多面函数代替传统的方法拟合海底地形的趋势项，通过算例验证了该方法的有效性，进一步提高了配置法在多波束测深数据处理中的推估精度。

图7 两种抗差方法测深点等价权值概况Fig.7 Equivalent weight of sounding points for two robust methods

(2) 针对多波束测深数据存在粗差或异常点时，常规最小二乘配置方法给出的协方差函数不能精确表征其统计特性的问题，提出了一种抗差最小二乘配置迭代解法，以迭代计算的方式给出了经验协方差函数的参数估计和最小二乘配置解，通过算例验证了本文方法在一定程度上能够克服多波束测深数据中粗差或异常点的影响，提高了推估的精度，相比于传统的抗差方法，本文抗差方法能够更为有效地识别出测深数据中的异常点。

最小二乘配置问题的关键是确定合理的协方差函数，接下来还需进一步研究协方差函数的确定方法，从而进一步提高最小二乘配置推估的精度。

致谢：感谢东华理工大学王胜平博士提供多波束测深数据。

[1] 朱庆, 李德仁. 多波束测深数据的误差分析与处理[J]. 武汉测绘科技大学学报, 1998, 23(1): 1-4, 46. ZHU Qing, LI Deren. Error Analysis and Processing of Multibeam Soundings[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1998, 23(1): 1-4, 46.

[2] 陆秀平, 黄谟涛, 翟国君, 等. 多波束测深数据处理关键技术研究进展与展望[J]. 海洋测绘, 2016, 36(4): 1-6, 11.LU Xiuping, HUANG Motao, ZHAI Guojun, et al. Development and Prospect of Key Technologies for Multibeam Echosounding Data Processing[J]. Hydrographic Surveying and Charting, 2016, 36(4): 1-6, 11.

[3] 王海栋, 柴洪洲, 王敏. 多波束测深数据的抗差Kriging拟合[J]. 测绘学报, 2011, 40(2): 238-242, 248. WANG Haidong, CHAI Hongzhou, WANG Min. Multibeam Bathymetry Fitting Based on Robust Kriging[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(2): 238-242, 248.

[4] BORRE K. Mathematical Foundation of Geodesy: Selected Papers of Torben Krarup[M]. Berlin Heidelberg: Springer, 2006.

[5] MORITZ H. Advanced Physical Geodesy[M]. Karlsruhe: Wichmann, Abacus Press, 1980.

[6] WOLF H.Über Verallgemeinerte Kollokation[J]. Zeitschrift für Vermessungswesen, 1974, 99(11): 475-478.

[7] EL-FIKY G S, KATO T. Continuous Distribution of the Horizontal Strain in the Tohoku District, Japan, Predicted by Least-squares Collocation[J]. Journal of Geodynamics, 1998, 27(2): 213-236.

[8] 柴洪洲, 崔岳, 明锋. 最小二乘配置方法确定中国大陆主要块体运动模型[J]. 测绘学报, 2009, 38(1): 61-65. DOI: 10.3321/j.issn:1001-1595.2009.01.011. CHAI Hongzhou,CUI Yue,MING Feng.The Determination of Chinese Mainland Crustal Movement Model Using Least-squares Collocation[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(1): 61-65. DOI: 10.3321/j.issn:1001-1595.2009.01.011.

[9] 江在森, 刘经南. 应用最小二乘配置建立地壳运动速度场与应变场的方法[J]. 地球物理学报, 2010, 53(5): 1109-1117. JIANG Zaisen, LIU Jingnan. The Method in Establishing Strain Field and Velocity Field of Crustal Movement Using Least Squares Collocation[J]. Chinese Journal of Geophysics, 2010, 53(5): 1109-1117.

[10] YANG Y X,ZENG A M,ZHANG J.Adaptive Collocation with Application in Height System Transformation[J]. Journal of Geodesy, 2009, 83(5): 403-410.

[11] JARMOLOWSKI W. Least Squares Collocation with Uncorrelated Heterogeneous Noise Estimated by Restricted Maximum Likelihood[J]. Journal of Geodesy, 2015, 89(6): 577-589.

[12] YOU R J, HWANG H W. Coordinate Transformation between Two Geodetic Datums of Taiwan by Least-Squares Collocation[J]. Journal of Surveying Engineering, 2006, 132(2): 64-70.

[13] 徐卫明, 刘雁春. 海道测量中非定位点的精密确定方法[J]. 系统仿真学报, 2007, 19(20): 4612-4615. XU Weiming, LIU Yanchun. Methods of Determining Non-Position Points in Hydrographic Survey[J]. Journal of System Simulation, 2007, 19(20): 4612-4615.

[14] 王海栋, 柴洪洲, 赵东明. 基于抗差最小二乘配置的海底地形生成研究[J]. 系统仿真学报, 2010, 22(9): 2091-2094, 2099. WANG Haidong, CHAI Hongzhou, ZHAO Dongming. Research on Seabed Terrain Generation Based on Robust Least-squares Collocation[J]. Journal of System Simulation, 2010, 22(9): 2091-2094, 2099.

[15] 王海栋. 多波束系统测深异常处理理论与方法研究[D]. 郑州: 信息工程大学, 2010. WANG Haidong.Research on the Theories and Algorithms of Processing the Bathymetry Outlier in Multibeam Echosounder System[D]. Zhengzhou: Information Engineering University, 2010.

[16] 薛树强, 杨元喜. 广义反距离加权空间推估法[J]. 武汉大学学报(信息科学版), 2013, 38(12): 1435-1439. XUE Shuqiang, YANG Yuanxi. Generalized Inverse Distance Weighting Method for Spatial Interpolation[J]. Geomatics and Information Science of Wuhan University, 2013, 38(12): 1435-1439.

[17] 黄维彬. 近代平差理论及其应用[M]. 北京: 解放军出版社, 1992: 92-118. HUANG Weibin. Theory of Adjustment and Its Applications in Geodesy[M]. Beijing: The PLA Press, 1992: 92-118.

[18] YANG Yuanxi. Robustfying Collocation[J]. Manuscripta Geodaetica, 1992, 17(1): 21-28.

[19] HARDY R L. Multiquadric Equations of Topography and Other Irregular Surfaces[J]. Journal of Geophysical Research, 1971, 76(8): 1905-1915.

[20] HARDY R L. The Application of Multiquadric Equations and Point Mass Anomaly Models to Crustal Movement Studies[R]. NOAA Technical Report NOS 76 NGS 11. Rockville, Maryland: US Department of Commerce, National Oceanic and Atmospheric Administration, National Ocean Survey, National Geodetic Survey, 1978.

[21] 谢曦霖, 许才军, 温扬茂, 等. 一种基于多面函数的改进最小二乘配置方法[J]. 武汉大学学报(信息科学版). DOI: 10.13203/j.whugis20150664. XIE Xilin, XU Caijun, WEN Yangmao, et al. A Refined Least Squares Collocation Method Based on Multiquadric Function[J]. Geomatics and Information Science of Wuhan University. DOI: 10.13203/j.whugis20150664.

[22] 杨元喜. 抗差估计理论及其应用[M]. 北京: 八一出版社, 1993. YANG Yuanxi. Robust Estimation Theory and Its Applications[M]. Beijing: Bayi Publishing House, 1993.

[23] SCHAFFRIN B, BOCK Y. Geodetic Deformation Analysis Based on Robust Inverse Theory[J]. Manuscripta Geodaetica, 1994, 19(1): 31-44.

[24] 刘念. 协方差函数的抗差拟合[J]. 测绘科学, 2001, 26(3): 25-28. LIU Nian. Robust Fitting of Covariance Function[J]. Science of Surveying and Mapping, 2001, 26(3): 25-28.

[25] BROUNS G, DE WULF A, CONSTALES D. Multibeam Data Processing: Adding and Deleting Vertices in a Delaunay Triangulation[J]. Hydrographic Journal, 2001(101): 3-9.

(责任编辑：陈品馨)

WANG Leyang (1983—), male, PhD, associate professor, majors in geodetic inversion and geodetic data processing.

Multi-beam Bathymetry Data Processing Using Iterative Algorithm of Robust Least Squares Collocation

WANG Leyang1,2,3,CHEN Hanqing1,2

1. Faculty of Geomatics, East China University of Technology, Nanchang 330013, China; 2. Key Laboratory of Watershed Ecology and Geographical Environment Monitoring, NASG, Nanchang 330013, China; 3. Key Laboratory for Digital Land and Resources of Jiangxi Province, Nanchang 330013, China

In the process of dealing with multi-beam bathymetry data by least squares collocation, the quadric curved mathematical model of trend term can not express accurately the whole variation trend of seafloor topography in general. Moreover, the covariance function estimated by general method is incapable of accurately expressing statistical characteristics with the multi-beam bathymetry data contains gross errors or outliers. So the iterative algorithm of robust least squares collocation is proposed in this paper. Firstly, the initial weight matrix of observations and the initial parameters of covariance function are both given in this method, then the trend term is fitted by polyhedral function and equivalent weights scheme is applied into robust estimation in this method. Finally, the robust parameters of covariance function and solutions of least squares collocation are iteratively calculated. The experimental results show that the method proposed in this paper can express well the whole variation trend of seafloor topography and overcome the effect of gross error or outlier in multi-beam bathymetry data to a certain extent. Compared with the conventional robust method, the proposed method in this paper more effectively probes the outliers in bathymetry data with the robust and better predicted results.

least squares collocation; seabed terrain generation; covariance function; robust; trend term

National Natural Science Foundation of China (Nos.41664001; 41204003); Support Program for Outstanding Youth Talents in Jiangxi Province (No.20162BCB23050); National Key Research and Development Program(No.2016YFB0501405); Science and Technology Project of the Education Department of Jiangxi Province (No.GJJ150595); the Project of Key Laboratory for Digital Land and Resources of Jiangxi Province(No.DLLJ201705).

王乐洋，陈汉清.多波束测深数据处理的抗差最小二乘配置迭代解法[J].测绘学报，2017,46(5)：658-665.

10.11947/j.AGCS.2017.20160491. WANG Leyang, CHEN Hanqing.Multi-beam Bathymetry Data Processing Using Iterative Algorithm of Robust Least Squares Collocation[J]. Acta Geodaetica et Cartographica Sinica,2017,46(5):658-665. DOI:10.11947/j.AGCS.2017.20160491.

2016-10-17

王乐洋(1983—)，男，博士，副教授，主要研究方向为大地测量反演及大地测量数据处理。

E-mail： wleyang@163.com

P229

1001-1595(2017)05-0658-08

国家自然科学基金 (41664001; 41204003)；江西省杰出青年人才资助计划项目(20162BCB23050)；国家重点研发计划(2016YFB0501405)；江西省教育厅科技项目(GJJ150595) ；江西省数字国土重点实验室开放研究基金(DLLJ201705)

修回日期： 2017-04-13