于家成　徐學(xué)強(qiáng)

【摘 要】提出了利用正態(tài)分布進(jìn)行多波束系統(tǒng)水深地圖重建的方法，首先根據(jù)二維正態(tài)分布對距離反比權(quán)重插值算法進(jìn)行改進(jìn)，形成正態(tài)插值算法，用該算法對原始水深數(shù)據(jù)進(jìn)行插值以模擬實(shí)際地形；繼而根據(jù)多波束系統(tǒng)測深原理對模擬地形進(jìn)行航向東偏北45°角的水深測量，通過5個(gè)校正因子對實(shí)時(shí)測量數(shù)據(jù)進(jìn)行校正；最后使用正態(tài)插值算法對改正后的測量數(shù)據(jù)重新插值，完成實(shí)時(shí)測量地圖的重建。使用某海域真實(shí)水深數(shù)據(jù)進(jìn)行仿真，結(jié)果表明，校正后的實(shí)時(shí)地圖相對直接測量誤差減小了3.3米，建圖精度獲得了明顯提高。

【關(guān)鍵詞】多波束數(shù)字測深；正態(tài)插值算法；地圖重建

Normal Interpolation Reconstruction of Multibeam Digital Bathymetry Map

YU Jia-cheng1 XU Xue-qiang2

（1. Department of Electronics and Information Engineering， North China Institute of Science and Technology ， Hebei Langfang， 065201， China；2. China South Industries Research Academy， Beijing，100089， China）

【Abstract】The method that realtime multibeam digital bathymetry map is reconstructed according to the law that seabed elevation follows normal distribution is proposed. At first， inverse distance weighted（IDW） is transformed to normal distribution weighted（NDW）， then the raw depth data is interpolated using NDW to imitate actual terrain. Secondly， the imitated terrain is measured along 45 degrees north east direction and realtime data are corrected by five factors； Thirdly， the corrected data are interpolated by NDW to construct realtime terrain map. The construction process is verified by using practical water depth data of certain seafloor. Simulation results show that the reconstructed map error is decreased by 3.3m relative to previous， thus the map reconstruction precision is improved obviously.

【Key words】Multibeam digital bathymetry；Normal distribution weighted（NDW）；Map reconstruction

多波束水深測量的關(guān)鍵在于建立高分辨率的海底數(shù)字地形模型（digital terrain model DTM）[1]， 這是目前多波束數(shù)據(jù)處理研究中的一個(gè)重點(diǎn)。為實(shí)現(xiàn)大數(shù)據(jù)量離散深度數(shù)據(jù)的快速海底數(shù)字地形建模和可視化，國內(nèi)外研究人員通常采用對離散深度數(shù)據(jù)進(jìn)行格網(wǎng)化處理[2]，其實(shí)質(zhì)是對深度數(shù)據(jù)進(jìn)行內(nèi)插，使離散分布的數(shù)據(jù)點(diǎn)變?yōu)橐?guī)則柵格分布。根據(jù)數(shù)據(jù)量的大小，一般分為中小規(guī)模數(shù)據(jù)的插值和大規(guī)模數(shù)據(jù)的插值大類[3]。本文以單平面換能器多波束系統(tǒng)為例建立仿真模型，對波束的立體寬度，從橫縱兩個(gè)方向?qū)λ顪y量的影響進(jìn)行改正，以提高水深測量精度，在此基礎(chǔ)上，采用以二維正態(tài)分布構(gòu)建權(quán)重的算法來構(gòu)建規(guī)則格網(wǎng)（GRID）的DTM。

1 距離反比權(quán)重法的改進(jìn)

1.1 算法原理

距離反比權(quán)重法（inverse distance weighted，IDW）基于如下的假定： 相互之間距離較近的對象相對于距離較遠(yuǎn)的對象之間有更強(qiáng)的相似性。在連續(xù)的自然地形表面上所采集的離散點(diǎn)高程屬性符合這種假定。IDW采用周邊已經(jīng)測定的采樣點(diǎn)（已知點(diǎn)）來確定未知點(diǎn)的數(shù)值， 距離該點(diǎn)位近的已知點(diǎn)對其數(shù)值的影響大于距離較遠(yuǎn)的點(diǎn)。因此，IDW假定數(shù)據(jù)集中的每一個(gè)已知點(diǎn)都對該未知點(diǎn)存在局部的影響力，這種影響力的大小隨著距離的增大而逐漸消失。IDW給距離較近的點(diǎn)分配一個(gè)比較大的權(quán)值，而給較遠(yuǎn)的點(diǎn)分配一個(gè)較小的權(quán)值以體現(xiàn)這種影響力[1]。即[4]

設(shè)空間待插點(diǎn)為點(diǎn)P（xP，yP，zP），P點(diǎn)鄰域內(nèi)有已知離散點(diǎn)Qi（xi，yi，zi），i=1，2，…，n，利用距離反比權(quán)重法對P點(diǎn)的屬性值z進(jìn)行插值。待插點(diǎn)的屬性值f（x，y）是待插點(diǎn)鄰域內(nèi)已知離散點(diǎn)屬性值的加權(quán)平均，

zp= w（d ）z / w（d ）（1）

式中，w（di）是權(quán)函數(shù)， 其值與待插點(diǎn)與鄰域內(nèi)離散點(diǎn)之間的距離有關(guān)， 是距離的k次方的倒數(shù)，w（di）=1/d2i。k值越低，插值結(jié)果越平滑，k值越高，插值結(jié)果細(xì)節(jié)越詳細(xì)，通常取為2。zi是第i個(gè)離散點(diǎn)的數(shù)值，di是該點(diǎn)到網(wǎng)格點(diǎn)之間的距離[5]。

1.2 算法改進(jìn)

大量研究表明：自然地形的起伏滿足隨機(jī)分布規(guī)律，而DEM是表示地面高程空間分布的數(shù)字地面模型， 反映自然地形的起伏變化，確實(shí)服從一定的統(tǒng)計(jì)分布規(guī)律，一般可以用正態(tài)分布代替[6]，即以未知點(diǎn)為中心，未知點(diǎn)與周圍點(diǎn)之間隨位置的遠(yuǎn)近呈現(xiàn)出一種正態(tài)分布的關(guān)系。未知點(diǎn)對自身的影響最大，具有最大的權(quán)值，周圍點(diǎn)隨著距離的增大權(quán)值逐漸減小，符合正態(tài)分布。對海底地形而言，點(diǎn)相對未知點(diǎn)所具有的權(quán)重為

w（d）=f（x，y）= exp = exp （2）

由于距離越遠(yuǎn)， 影響力越小。隨著距離的增大，己知點(diǎn)對于未知點(diǎn)的影響實(shí)際上已經(jīng)可以忽略不計(jì)， 為了加快運(yùn)算速度， 設(shè)定這些點(diǎn)的權(quán)重為0。因此， 需要指定一個(gè)未知點(diǎn)的搜索鄰域， 只考慮位于該鄰域內(nèi)的已知點(diǎn)對未知點(diǎn)所造成的影響。鄰域的形狀與數(shù)據(jù)集的特性有關(guān)： 如果權(quán)值確定不需要考慮方向性的影響， 則鄰域的形狀可以是一個(gè)圓形[6]

1.3 搜索半徑的確定

在文獻(xiàn)[7]中，使用一矩形作為搜索區(qū)域，對落在矩形窗內(nèi)的有效數(shù)據(jù)，都進(jìn)行插值。為了防止失真，設(shè)立有效數(shù)據(jù)個(gè)數(shù)門限，當(dāng)矩形窗內(nèi)的有效數(shù)據(jù)個(gè)數(shù)小于此下限值，則此網(wǎng)格節(jié)點(diǎn)為無效值。這個(gè)門限缺省值為1，即只有矩形窗內(nèi)存在數(shù)據(jù)，該網(wǎng)格節(jié)點(diǎn)即為有效點(diǎn)，但是為了確保插值的準(zhǔn)確性，推薦值為3。在本文中，為了符合正態(tài)權(quán)重公式w（d）= exp ，采用一圓形區(qū)域，如果所用網(wǎng)格為不完全規(guī)則的矩形網(wǎng)格，為了兼顧運(yùn)算速度和確保插值的準(zhǔn)確，把半徑取為最大網(wǎng)格的對角線長度，與文獻(xiàn)[1]不同的是，網(wǎng)格點(diǎn)為已知數(shù)據(jù)點(diǎn)，格內(nèi)的任一點(diǎn)為待插點(diǎn)，這樣任意待插點(diǎn)都至少有4個(gè)已知點(diǎn)為其插值作出貢獻(xiàn)，也為網(wǎng)格大小不一時(shí)保留一定的余度，以防在待插點(diǎn)周圍搜索不到已知點(diǎn)，如圖1所示，|AB|=max（經(jīng)度網(wǎng)格間距），|AC|=max（緯度網(wǎng)格間距），取搜索半徑R=|AD|。如果最大經(jīng)度網(wǎng)格間距與最大緯度網(wǎng)格間距不從屬同一個(gè)網(wǎng)格，仍然這樣取搜索半徑。

2 多波束系統(tǒng)水深掃描

2.1 多波束數(shù)字測深模型

如圖2所示，設(shè)沿船體橫向（垂直龍骨方向）第i個(gè)波束投射到海底到換能器所處方向的橫向傾角為β，沿船體運(yùn)動方向（平行龍骨方向）的縱向傾角為，換能器所在位置P。設(shè)船體從原點(diǎn)開始沿方向 航行，航速V，采樣時(shí)間T，由于存在橫縱兩個(gè)方向的波束寬度，第i個(gè)波束照射的海底離散點(diǎn)形成一個(gè)矩形點(diǎn)集Vi。

2.2 測深搜索算法設(shè)計(jì)

設(shè)船體沿以角度東偏北45°，船速6m/s，勻速直線航行50s，采樣時(shí)間1s。在進(jìn)行地圖掃描時(shí)，首先把換能器中心歸結(jié)到最近網(wǎng)格點(diǎn)，以該點(diǎn)為中心，以2max|z|·tg（iθ）為邊長建立正方形作為搜索區(qū)域，搜索點(diǎn)數(shù)為int +1，（grid為網(wǎng)格間距）；得到ri（xn，ym），從而得到Zi、Li，腳印位置（xi，yi），按照 方向，以足夠大矩形區(qū)域進(jìn)行搜索，中間波束腳印位置進(jìn)行最近網(wǎng)格點(diǎn)歸結(jié)后作為向外邊緣的下一個(gè)波束的搜索區(qū)域的起點(diǎn)。搜索完成后產(chǎn)生的16×51水深數(shù)據(jù)地圖如圖3所示，由該圖可以發(fā)現(xiàn)，測船掃描過一段左低右高的慢坡曲面。水深在-980米～-1060米之間。

3 水深測量校正

多波束系統(tǒng)的測深換能器發(fā)射的聲波具有一個(gè)立體波束角，聲波到達(dá)海底或碰到不同的介質(zhì)界面時(shí)返射，再被換能器接收，接收換能器收到的信號可能是在這個(gè)波束角所覆蓋的任何水深。正確深度應(yīng)該是波束角所覆蓋區(qū)域的中心點(diǎn)水深值，顯然，它不同于聲波覆蓋區(qū)域的邊緣點(diǎn)水深值，因此產(chǎn)生水深測量誤差。此項(xiàng)誤差在海底水深變化較大的地方尤為明顯[4]，尤其在實(shí)測水深較深或海底坡度較大的地形時(shí)，由于存在波束寬度，往往產(chǎn)生地形圖等深線向斷面方向前后凸凹、彎曲的情形[8]，形成虛假地形，從而使繪制的海底地圖與實(shí)際地形存在差異。在水深測量中，是以最短距離作為記錄長度，由于存在波束寬度及地形的起伏，并且隨著波束投射角的增大，測量的誤差越來越大，尤其邊緣波束更趨明顯。在以往的實(shí)際工作中，一般認(rèn)為約占10%的外側(cè)波束（左右各占5%）定義為邊緣波束，為了保證數(shù)據(jù)的精度滿足要求，常常強(qiáng)制性地將這些波束舍棄[9]。為了減小測量誤差，需要對測得的數(shù)據(jù)進(jìn)行地形補(bǔ)償。

在多波束水深測量中，由于存在橫向波束寬度，隨著如圖2所示中α，β的影響。為了減小測量誤差，需要對測得的數(shù)據(jù)進(jìn)行地形補(bǔ)償。以因子A1=cos（（i-1）θ-α）·sec（（i- ）θ-α），A2=cos（α-iθ）·sec（α-（i- ）θ），A3=sec（|（i- ）θ-α|）來校正沿船體橫向（垂直龍骨）的傾斜角；以B1= ，B2= 來校正沿船體縱向（平行龍骨）的傾斜角。校正前后的水深誤差如圖4（a）（b）所示，比較發(fā)現(xiàn)，校正后的水深誤差極大減小，精度提高4倍左右。

4 實(shí)時(shí)圖重建

4.1 實(shí)時(shí)圖重插值

對經(jīng)過校正后的水深值，按照校正后的水深進(jìn)行插值，確定插值區(qū)域的原則是，插值區(qū)域?yàn)閷?shí)時(shí)掃描圖的中間矩形區(qū)域，并且最大程度的覆蓋掃描區(qū)。在掃描過程中，相鄰波束腳印之間的距離往往較大，并且隨著水深變化、地形起伏也不斷變化，在建立實(shí)時(shí)掃描地圖時(shí)需要在中間進(jìn)行插值填補(bǔ)以連接起來，沿航行方向及其垂直方向建立規(guī)則格網(wǎng)，如圖5中所示矩形區(qū)域，從換能器所在航線的中點(diǎn)開始，向左右、上下蔓延，為保持一致，每次蔓延的網(wǎng)格間距與初始插值時(shí)相同（實(shí)際插值時(shí)沒此限制，可以根據(jù)需要設(shè)置）。使用ceil取整，以使得待插點(diǎn)都是整數(shù)個(gè)網(wǎng)格間距，同時(shí)又在插值區(qū)域內(nèi)，其中，“×”為已知腳印點(diǎn)，格網(wǎng)的交叉點(diǎn)為待插值點(diǎn)。

使用本文1中的方法計(jì)算校正后水深的平均值和標(biāo)準(zhǔn)差，進(jìn)行正態(tài)插值，建立新的水深實(shí)時(shí)測量圖如圖6（a）所示，相應(yīng)的建圖誤差如圖6（b）所示。

4.2 重建誤差分析

將地圖重建后的水深值、水深測量值、測量校正值，分別與同位置原始水深進(jìn)行比較，同位置原始水深值通過對初始10×10格網(wǎng)進(jìn)行正態(tài)插值獲得，它們之間的誤差如表1所示，由表1發(fā)現(xiàn)，改正后的平均誤差與最大誤差都顯著減小，地圖重建后的平均誤差及最大誤差相比測量改正后的誤差進(jìn)一步減小，這主要是因?yàn)樵谶x取插值區(qū)域時(shí)選擇了掃描地區(qū)的中間矩形區(qū)域，邊緣波束被舍掉，而邊緣波束測量及改正后的誤差都在整個(gè)掃描圖中最大所致。

5 結(jié)論

本文通過正態(tài)權(quán)重的距離反比權(quán)重插值算法，對某海域的已有實(shí)際水深數(shù)據(jù)進(jìn)行插值，繼而使用多波束進(jìn)行掃描、校正，完成實(shí)時(shí)地圖的重建。海底地形變化多樣，細(xì)節(jié)復(fù)雜，使用5個(gè)因子的校正是近似的。同時(shí)，為了航行安全起見，可以保留適度的淺誤差，但如果取得的實(shí)時(shí)圖是為了匹配，則要提高精度，要根據(jù)實(shí)際情況，二者兼顧。

【參考文獻(xiàn)】

[1]馬建林，金菁，來向華.多波束測深海底數(shù)字地形模型的建立.海洋測繪，2005，25（5）：15-17.

[2]陸丹，李海森，魏玉闊，鄧平.多波束測深系統(tǒng)中的海底地形可視化技術(shù)研究[J].儀器儀表學(xué)報(bào)，2012，33（2）：450-456.

[3]何磊，袁德奎，陶建華.一種適用波浪潮流計(jì)算的地形插值算法[J].港工技術(shù)，2004，（4）：4-6.

[4]李章林，張夏林.距離平方反比法礦產(chǎn)資源儲量計(jì)算模塊設(shè)計(jì)與實(shí)現(xiàn)[J].地質(zhì)與勘探，2007，43（6）：92-97.

[5]靳國棟，劉衍聰，牛文杰.距離加權(quán)反比插值法和克里金插值法的比較[J].長春工業(yè)大學(xué)學(xué)報(bào)，2003，24（3）：53-57.

[6]劉曉萌，常占強(qiáng).數(shù)字高程模型數(shù)據(jù)分布規(guī)律的研究[J].河北師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2006，30（4）：473-477.

[7]高佩蘭，關(guān)永賢.多波束海底地形圖的制作方法探討[J].海洋測繪，2012，32（4）：71-73.

[8]褚宏憲，周小明，史慧杰，馮京.水深測量誤差分析與改正[J].物探與化探，2011，35（3）：358-362.

[9]關(guān)永賢，劉勝旋，劉方蘭.多波束測深系統(tǒng)邊緣波束的處理技術(shù)研究[J].中國人口·資源與環(huán)境，2003，13（68）：6-7.

[責(zé)任編輯：朱麗娜]