朱振華

(廣州海洋地質(zhì)調(diào)查局，廣州 510760)

1 概述

剔除了明顯粗差的河道多波束測量數(shù)據(jù)中，其邊緣波束尤其近岸測量條帶邊緣波束由于缺乏校正手段，受到的系統(tǒng)誤差影響較大。為了提高該部分?jǐn)?shù)據(jù)精度，必須對其系統(tǒng)誤差和小粗差進(jìn)行處理。由于系統(tǒng)誤差的影響，本來的正確數(shù)據(jù)可能被認(rèn)為是小粗差，原本是小粗差的點(diǎn)可能被遮蔽而作為正確數(shù)據(jù)，因此，在數(shù)據(jù)存在若干人工檢測點(diǎn)的情況下，本文采用半(非)參數(shù)法先處理系統(tǒng)誤差并修正高程值后再以抗差M估計(jì)法檢測異常，期望在最大程度上還原真實(shí)水深值。

2 半(非)參數(shù)模型

在內(nèi)河聲納測量中，系統(tǒng)誤差的主要來源是測船橫搖。由于聲納測量屬于動態(tài)測量，觀測結(jié)果存在著無法參數(shù)化的系統(tǒng)誤差。半?yún)?shù)模型為解決此類問題提供了可行之道，該模型表達(dá)為間接平差模型的延伸:

式中:L和Δ為n維觀測向量，Δ～(0，σ2)是觀測誤差向量，x為參數(shù)向量，B為列滿秩矩陣，S=[s1，s2，…，sn]T是一非參數(shù)待定量，作為系統(tǒng)誤差的描述。列出其誤差方程為:

由最小二乘原理VTPV=min求得的法方程只有t個方程，無法確定S和Δ的估值。根據(jù)Fisher和Hegland(1999)的觀點(diǎn)[1-2]，給出補(bǔ)償最小二乘準(zhǔn)則為:

式中:R為正規(guī)化矩陣，α稱為規(guī)范化參數(shù)或平滑因子。以式(2)為條件，由拉格朗日乘常數(shù)法構(gòu)造函數(shù):令，經(jīng)過一系列推導(dǎo)可得:

其中 N=BTPB，M=P+αR-PBN-1BTP。由文獻(xiàn)[1]和[2]的觀點(diǎn)，取:

G的定義請見上述文獻(xiàn)。由于rank(R=n-1＜n)，需要增加一個約束條件再次由拉格朗日乘常數(shù)法構(gòu)造函數(shù)求條件極值。當(dāng)測船在變化平緩的河道以較慢速度測量時，可以認(rèn)為一個橫搖期間所測量的河床點(diǎn)實(shí)際高程相等，而橫搖是周期性運(yùn)動，在一次橫搖中對高程值的正負(fù)影響基本相互抵消，在一定范圍測區(qū)內(nèi)可將系統(tǒng)誤差之和看作零，即，得:

拉格朗日函數(shù)為:

令W=P-PBN-1BTP，推導(dǎo)得到式(8)的解:

確定S的估值和正規(guī)化矩陣R后，采用L-曲線法確定α的值。以一小段河道為例，當(dāng)在近岸處人工測量若干深度并作為真值)z代替式(1)中的Bx，即以檢測點(diǎn)的真誤差l=z-)z為觀測值時，系統(tǒng)誤差為S，觀測噪聲為Δ，則可以將式(1)的系數(shù)陣看作零，這時稱為非參數(shù)估計(jì)，式(1)可表達(dá)為:

分別由式(10)、式(11)可推出:

其中 M=PαR，γ =eTM-1WL/(eTM-1e)，W=P=E。對由式(13)確定的系統(tǒng)誤差向量S，建立如下函數(shù)關(guān)系:

將該河段邊緣波束共數(shù)千個原始數(shù)據(jù)代入上式，即可求得各點(diǎn)的系統(tǒng)誤差。表1中x為橫坐標(biāo)，y為縱坐標(biāo)，為方便起見，縱橫坐標(biāo)分別省去前面三、四位數(shù)字，即各點(diǎn)原坐標(biāo)應(yīng)為(455***.**，3799***，**)，z為多波束測深值，)z為同一平面坐標(biāo)的人工測量值，共21個。S1為非參數(shù)法去除的系統(tǒng)誤差。

3 抗差M估計(jì)濾波法

根據(jù)黃謨濤的觀點(diǎn)[3]，基于M估計(jì)的抗差解為:

表1 兩種方法去除的系統(tǒng)誤差對比 m

第k+1步迭代解為:

式中:A為方程系數(shù)矩陣，L是自由項(xiàng)，X為模型待定參數(shù)向量，pi為觀測值的權(quán)，為等價權(quán)陣。等價權(quán)值pi采用楊元喜提出的 IGGⅢ[4]方案:

式中:v′i=vi/，k0和 k1分別取 1.0 ～1.5 和 2.0 ～3.0，權(quán)值由上一次迭代得到的加權(quán)平均值計(jì)算殘差，再判斷比較條件來確定。k次迭代后的為:

由迭代解表達(dá)式，設(shè)受檢測點(diǎn)鄰域內(nèi)某點(diǎn)深度為zi(i=1，…，m)，可得基于抗差M估計(jì)的選權(quán)迭代加權(quán)平均值計(jì)算公式如下:

4 結(jié)論與建議

去除系統(tǒng)誤差后，高程的歪曲現(xiàn)象得以暴露。高程的歪曲包括粗差遮蔽和高程放大，即去除系統(tǒng)誤差前，某些點(diǎn)會被錯誤判斷為異?；蛟緦儆趹岩僧惓｜c(diǎn)卻被漏檢。系統(tǒng)誤差被去除后，點(diǎn)的真值最近似值就得以還原。因此，沿用2節(jié)中的例子，對去除系統(tǒng)誤差后的示例數(shù)據(jù)進(jìn)行小異常的檢測。取所有迭代權(quán)穩(wěn)定后為迭代停止條件，用抗差M估計(jì)法篩選異常值。表2是去除系統(tǒng)誤差前后提取異常值的對比，而圖1是去除系統(tǒng)誤差前后實(shí)驗(yàn)河段原始數(shù)據(jù)及異常值分布圖，借助 MATLAB 顯示[6]。

表2中可以看到，去除系統(tǒng)誤差前篩選到8個異常值，去除系統(tǒng)誤差后篩選到9個異常值，在去除系統(tǒng)誤差前后篩選的所有異常值中，前7個異常值都被檢測到，它們的平面坐標(biāo)一致。去除系統(tǒng)誤差前篩選到的第8個點(diǎn)(加粗顯示)因受到系統(tǒng)誤差影響導(dǎo)致高程被放大而被誤列為異常點(diǎn)，該點(diǎn)去除系統(tǒng)誤差后高程為30.05 m;去除系統(tǒng)誤差后篩選到的第8、9個異常點(diǎn)(加粗顯示)在去除系統(tǒng)誤差前因受到系統(tǒng)誤差遮蔽而未被檢測到，其實(shí)際高程值分別為30.17 m和30.39 m。圖1中(a)和(b)中黑色點(diǎn)分別為表2中所篩選的異常值。

從上述結(jié)果可以看出，本文所采用的兩種方法在處理河道多波束測量數(shù)據(jù)高程誤差歪曲現(xiàn)象中取得了良好的效果，同時有效地修正了高程值。但是該處理方法也存在一定的缺點(diǎn)，目前只在變化平緩的河道可以取得較好的效果，應(yīng)進(jìn)一步改善算法，使其局限性得以突破。

表2 去除系統(tǒng)誤差前后篩選的異常值

圖1 去除系統(tǒng)誤差前后剔除的異常值示意圖

[1] 胡昌宏.半?yún)?shù)模型的估計(jì)方法及其應(yīng)用[D].武漢大學(xué)博士學(xué)位論文，2004.

[2] 孫海燕，吳云.半?yún)?shù)回歸與模型精化[J].武漢大學(xué)學(xué)報(bào)信息科學(xué)版.2002，27(2):172-174.

[3] 黃謨濤，翟國君，王瑞，歐陽永忠，管錚.海洋測量中異常數(shù)據(jù)的定位研究[J]. 海洋測繪，1999，(2):10-19.

[4] 楊元喜.抗差估計(jì)理論及其應(yīng)用[M].北京:八一出版社，1993.

[5] 趙建虎.多波束深度及圖像數(shù)據(jù)處理方法研究[D].武漢大學(xué)博士學(xué)位論文，2000.

[6] 王紅，邸國輝，熊燕.水下地形測量高程異常點(diǎn)剔除方法研究[J]. 四川測繪，2004，27(1):36-38.