朱振華
(廣州海洋地質(zhì)調(diào)查局,廣州 510760)
剔除了明顯粗差的河道多波束測量數(shù)據(jù)中,其邊緣波束尤其近岸測量條帶邊緣波束由于缺乏校正手段,受到的系統(tǒng)誤差影響較大。為了提高該部分數(shù)據(jù)精度,必須對其系統(tǒng)誤差和小粗差進行處理。由于系統(tǒng)誤差的影響,本來的正確數(shù)據(jù)可能被認為是小粗差,原本是小粗差的點可能被遮蔽而作為正確數(shù)據(jù),因此,在數(shù)據(jù)存在若干人工檢測點的情況下,本文采用半(非)參數(shù)法先處理系統(tǒng)誤差并修正高程值后再以抗差M估計法檢測異常,期望在最大程度上還原真實水深值。
在內(nèi)河聲納測量中,系統(tǒng)誤差的主要來源是測船橫搖。由于聲納測量屬于動態(tài)測量,觀測結(jié)果存在著無法參數(shù)化的系統(tǒng)誤差。半?yún)?shù)模型為解決此類問題提供了可行之道,該模型表達為間接平差模型的延伸:
式中:L和Δ為n維觀測向量,Δ~(0,σ2)是觀測誤差向量,x為參數(shù)向量,B為列滿秩矩陣,S=[s1,s2,…,sn]T是一非參數(shù)待定量,作為系統(tǒng)誤差的描述。列出其誤差方程為:
由最小二乘原理VTPV=min求得的法方程只有t個方程,無法確定S和Δ的估值。根據(jù)Fisher和Hegland(1999)的觀點[1-2],給出補償最小二乘準則為:
式中:R為正規(guī)化矩陣,α稱為規(guī)范化參數(shù)或平滑因子。以式(2)為條件,由拉格朗日乘常數(shù)法構(gòu)造函數(shù):令,經(jīng)過一系列推導(dǎo)可得:
其中 N=BTPB,M=P+αR-PBN-1BTP。由文獻[1]和[2]的觀點,取:
G的定義請見上述文獻。由于rank(R=n-1<n),需要增加一個約束條件再次由拉格朗日乘常數(shù)法構(gòu)造函數(shù)求條件極值。當(dāng)測船在變化平緩的河道以較慢速度測量時,可以認為一個橫搖期間所測量的河床點實際高程相等,而橫搖是周期性運動,在一次橫搖中對高程值的正負影響基本相互抵消,在一定范圍測區(qū)內(nèi)可將系統(tǒng)誤差之和看作零,即,得:
拉格朗日函數(shù)為:
令W=P-PBN-1BTP,推導(dǎo)得到式(8)的解:
確定S的估值和正規(guī)化矩陣R后,采用L-曲線法確定α的值。以一小段河道為例,當(dāng)在近岸處人工測量若干深度并作為真值)z代替式(1)中的Bx,即以檢測點的真誤差l=z-)z為觀測值時,系統(tǒng)誤差為S,觀測噪聲為Δ,則可以將式(1)的系數(shù)陣看作零,這時稱為非參數(shù)估計,式(1)可表達為:
分別由式(10)、式(11)可推出:
其中 M=PαR,γ =eTM-1WL/(eTM-1e),W=P=E。對由式(13)確定的系統(tǒng)誤差向量S,建立如下函數(shù)關(guān)系:
將該河段邊緣波束共數(shù)千個原始數(shù)據(jù)代入上式,即可求得各點的系統(tǒng)誤差。表1中x為橫坐標,y為縱坐標,為方便起見,縱橫坐標分別省去前面三、四位數(shù)字,即各點原坐標應(yīng)為(455***.**,3799***,**),z為多波束測深值,)z為同一平面坐標的人工測量值,共21個。S1為非參數(shù)法去除的系統(tǒng)誤差。
根據(jù)黃謨濤的觀點[3],基于M估計的抗差解為:
表1 兩種方法去除的系統(tǒng)誤差對比 m
第k+1步迭代解為:
式中:A為方程系數(shù)矩陣,L是自由項,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)平均值計算殘差,再判斷比較條件來確定。k次迭代后的為:
由迭代解表達式,設(shè)受檢測點鄰域內(nèi)某點深度為zi(i=1,…,m),可得基于抗差M估計的選權(quán)迭代加權(quán)平均值計算公式如下:
去除系統(tǒng)誤差后,高程的歪曲現(xiàn)象得以暴露。高程的歪曲包括粗差遮蔽和高程放大,即去除系統(tǒng)誤差前,某些點會被錯誤判斷為異?;蛟緦儆趹岩僧惓|c卻被漏檢。系統(tǒng)誤差被去除后,點的真值最近似值就得以還原。因此,沿用2節(jié)中的例子,對去除系統(tǒng)誤差后的示例數(shù)據(jù)進行小異常的檢測。取所有迭代權(quán)穩(wěn)定后為迭代停止條件,用抗差M估計法篩選異常值。表2是去除系統(tǒng)誤差前后提取異常值的對比,而圖1是去除系統(tǒng)誤差前后實驗河段原始數(shù)據(jù)及異常值分布圖,借助 MATLAB 顯示[6]。
表2中可以看到,去除系統(tǒng)誤差前篩選到8個異常值,去除系統(tǒng)誤差后篩選到9個異常值,在去除系統(tǒng)誤差前后篩選的所有異常值中,前7個異常值都被檢測到,它們的平面坐標一致。去除系統(tǒng)誤差前篩選到的第8個點(加粗顯示)因受到系統(tǒng)誤差影響導(dǎo)致高程被放大而被誤列為異常點,該點去除系統(tǒng)誤差后高程為30.05 m;去除系統(tǒng)誤差后篩選到的第8、9個異常點(加粗顯示)在去除系統(tǒng)誤差前因受到系統(tǒng)誤差遮蔽而未被檢測到,其實際高程值分別為30.17 m和30.39 m。圖1中(a)和(b)中黑色點分別為表2中所篩選的異常值。
從上述結(jié)果可以看出,本文所采用的兩種方法在處理河道多波束測量數(shù)據(jù)高程誤差歪曲現(xiàn)象中取得了良好的效果,同時有效地修正了高程值。但是該處理方法也存在一定的缺點,目前只在變化平緩的河道可以取得較好的效果,應(yīng)進一步改善算法,使其局限性得以突破。
表2 去除系統(tǒng)誤差前后篩選的異常值
圖1 去除系統(tǒng)誤差前后剔除的異常值示意圖
[1] 胡昌宏.半?yún)?shù)模型的估計方法及其應(yīng)用[D].武漢大學(xué)博士學(xué)位論文,2004.
[2] 孫海燕,吳云.半?yún)?shù)回歸與模型精化[J].武漢大學(xué)學(xué)報信息科學(xué)版.2002,27(2):172-174.
[3] 黃謨濤,翟國君,王瑞,歐陽永忠,管錚.海洋測量中異常數(shù)據(jù)的定位研究[J]. 海洋測繪,1999,(2):10-19.
[4] 楊元喜.抗差估計理論及其應(yīng)用[M].北京:八一出版社,1993.
[5] 趙建虎.多波束深度及圖像數(shù)據(jù)處理方法研究[D].武漢大學(xué)博士學(xué)位論文,2000.
[6] 王紅,邸國輝,熊燕.水下地形測量高程異常點剔除方法研究[J]. 四川測繪,2004,27(1):36-38.