基于無人船系統(tǒng)的水下DEM插值模型分析

秦世偉, 孫 炎, 武亞軍

(上海大學(xué)土木工程系, 上海 200444)

數(shù)字高程模型(digital elevation model, DEM)為一定區(qū)域內(nèi)連續(xù)變化的地表高程的數(shù)學(xué)表示, 是地理信息系統(tǒng)中最為重要的數(shù)據(jù)系統(tǒng)[1]．內(nèi)河航運(yùn)作為我國交通運(yùn)輸體系的重要組成部分,其水下DEM能有效反映地形的基本空間分布狀態(tài), 是航道正常運(yùn)轉(zhuǎn)的重要指標(biāo)之一．河道DEM豐富的地理信息,常被用于水土保持、水文分析、淤泥清理、環(huán)境整治等方面[2]．獲取河道DEM的傳統(tǒng)方法是利用六分儀、經(jīng)緯儀、測深桿等工具對(duì)水下地形點(diǎn)進(jìn)行測量, 這種方法不僅效率低、精度差,對(duì)作業(yè)人員的安全也會(huì)構(gòu)成一定的威脅．

近年來, 無人船(unmanned surface vessel, USV)系統(tǒng)在水下地形測繪中得到了廣泛應(yīng)用,它不僅能有效提高測量數(shù)據(jù)的精度,還能實(shí)現(xiàn)水域高程的全覆蓋測繪．但無法獲取河道DEM中所有地形點(diǎn)的信息,在實(shí)際應(yīng)用中,須利用插值模型估算未知地形點(diǎn)高程, 對(duì)DEM進(jìn)行二次加密;因此,選取合適的插值模型很大程度上影響著DEM的精確性．插值模型適用性的研究受到國內(nèi)外學(xué)者的關(guān)注．Arun[3]使用五種插值方法對(duì)印度博帕爾及其周邊地區(qū)的DEM進(jìn)行插值處理,發(fā)現(xiàn)在大多數(shù)情況下,克里金法的效果優(yōu)于其他插值方法; Nistor等[4]依據(jù)現(xiàn)有DEM數(shù)據(jù),采用反距離權(quán)重法和普通克里金法預(yù)測了新加坡四個(gè)地層的地下水位分布狀況; Diaconu等[5]對(duì)比反距離權(quán)重法和簡單克里金法對(duì)多瑙河下游觀測數(shù)據(jù)的插值結(jié)果,發(fā)現(xiàn)兩種模型的插值誤差較為接近; 徐洋等[6]采用六種插值模型對(duì)萊州灣和煙臺(tái)港的水下岸坡地形進(jìn)行插值,發(fā)現(xiàn)自然領(lǐng)域法和ANUDEM在不同尺度下的插值效果較好．上述研究僅關(guān)注了插值模型構(gòu)建DEM的實(shí)際應(yīng)用和精度對(duì)比,忽略了地形空間分布對(duì)插值模型精度的影響．

本文基于無人船系統(tǒng)在上海潘涇地區(qū)的實(shí)測數(shù)據(jù), 擬采用反距離權(quán)重法(inverse distance weighting, IDW)、普通克里金法(ordinary Kriging, OK)、移動(dòng)最小二乘法(moving least squares, MLS)、徑向基函數(shù)插值(radial basis function, RBF)、不規(guī)則三角網(wǎng)插值(triangulated irregular network, TIN)等五種插值模型構(gòu)建河道DEM,通過交叉驗(yàn)證和模型對(duì)比等方法,評(píng)價(jià)不同水下地形的插值模型精度,得出測區(qū)不同地表粗糙度的最優(yōu)插值模型,定性分析不同插值模型構(gòu)建水下DEM的特點(diǎn),并在此基礎(chǔ)上,提出一種優(yōu)化插值流程的改進(jìn)算法．

1 研究方法

1.1 研究區(qū)域與數(shù)據(jù)采集

潘涇是位于上海市潘涇路東側(cè)的一條自然河流, 南起荻涇, 向北流經(jīng)羅南、羅店等地區(qū),全長約19 km, 是寶山區(qū)重要的內(nèi)河航道之一．本文使用智能無人船系統(tǒng)(安徽科微智能科技有限公司, 合肥)搭載水平精度為10 mm±1 ppm, 垂直精度為20 mm±1 ppm的實(shí)時(shí)動(dòng)態(tài)測量系統(tǒng)(real-time kinematic, RTK)(上海華測導(dǎo)航技術(shù)股份有限公司, 上海)與精度為深度×0.1%±1 cm的單波束聲學(xué)測深儀(南京元厚電氣有限公司, 南京)在潘涇地區(qū)進(jìn)行水下地形點(diǎn)實(shí)測．本研究所用無人船及實(shí)測狀況如圖1所示．

圖1 無人船在潘涇實(shí)地測量Fig.1 USV used to obtain the underwater elevation of the Panjing River

采用雙次測量法驗(yàn)證無人船單波束測量精度,首次測量的斷面線設(shè)置為主測線, 二次測量的斷面線作為檢測線,將檢測線與主測線的交叉點(diǎn)水深值進(jìn)行比對(duì),共獲得交叉驗(yàn)證點(diǎn)4 132個(gè), 對(duì)照《水運(yùn)工程測量規(guī)范》(JTS 131—2012)的相關(guān)要求,在水深小于20 m的區(qū)域內(nèi), 深度誤差應(yīng)小于0.2 m, 分析發(fā)現(xiàn)共有37處超出允許誤差, 點(diǎn)位合格率為99.18%,精度符合要求．本文以主測線測量數(shù)據(jù)作為研究對(duì)象,剔除數(shù)據(jù)中的RTK非固定解,并根據(jù)中值濾波法自動(dòng)修正水深數(shù)據(jù)的突變異常點(diǎn),最終得到測區(qū)整體地形描述參數(shù)與實(shí)測高程點(diǎn)分布情況,如表1所示．

表1 地形描述參數(shù)與實(shí)測高程點(diǎn)分布Tab.1 Topographic parameters and the topographic distribution of the measuring points

1.2 插值模型

OK與IDW類似,也是通過區(qū)域內(nèi)實(shí)測點(diǎn)加權(quán)求和來估計(jì)插值點(diǎn)的高程,但不同于IDW的是, OK插值前須考慮實(shí)測點(diǎn)在空間位置上空間屬性的變異分布來求取變異函數(shù),變異函數(shù)須同時(shí)滿足無偏性和估計(jì)值最小的原則．常見的變異函數(shù)有高斯模型、指數(shù)模型、球狀模型等,運(yùn)算時(shí)首先要選取合適的變異函數(shù),然后將區(qū)域內(nèi)實(shí)測點(diǎn)高程進(jìn)行加權(quán)即可求取插值點(diǎn)的高程值．

圖2 RBF工作原理Fig.2 Algorithmic principle of RBF

RBF由輸入層、隱含層、輸出層三層神經(jīng)元構(gòu)成, 其結(jié)構(gòu)如圖2所示．RBF能夠在確保插值曲面通過所有點(diǎn)的同時(shí)保持曲率最小, 是一系列精確插值方法的組合[8],其輸入層和輸出層都是線性層,而隱含層使用徑向基核函數(shù)．常用的核函數(shù)有高斯核函數(shù)、多重二次曲面函數(shù)、薄板樣條曲面函數(shù)、自然三次樣條曲面函數(shù)等,核函數(shù)直接影響插值精度．

TIN通過將區(qū)域內(nèi)一系列離散的高程點(diǎn)用互不連接、互不重疊的三角形連接, 用連續(xù)的三角形近似模擬三維地形．TIN可在保持原有數(shù)據(jù)精度的同時(shí),通過不同層次的分辨率來逼近真實(shí)地形．生成三角網(wǎng)后,遍歷所有三角形柵格以判斷插值點(diǎn)所在的三角形區(qū)域, 利用重心坐標(biāo)插值即可估算得到插值點(diǎn)的高程．

1.3 評(píng)價(jià)方法

2 結(jié)果與分析

2.1 模型精度評(píng)價(jià)

為了保證插值精度的可比性, 依據(jù)測線布設(shè)情況將測區(qū)每隔50 m斷面劃分為一個(gè)子區(qū)段, 選取其中采樣點(diǎn)密度分布均勻的40個(gè)樣本區(qū)段, 利用ArcGIS地統(tǒng)計(jì)模塊中的Create Subsets功能,自動(dòng)抽取85%的樣本點(diǎn)作為高程內(nèi)插樣本數(shù)據(jù)集合,另外15%的樣本點(diǎn)作為驗(yàn)證數(shù)據(jù)集, 利用交叉檢驗(yàn)評(píng)價(jià)模型的插值精度．根據(jù)內(nèi)河河道起伏、切割程度不同的特點(diǎn),選取地表粗糙度Sr作為區(qū)段地形評(píng)價(jià)指標(biāo), 能夠在宏觀上體現(xiàn)地形起伏和侵蝕程度的變化[9], 計(jì)算公式為Sr=1/cosJ, 其中J為區(qū)段內(nèi)地形的平均坡度, 通過ArcGIS的坡度工具計(jì)算得到．在得到各區(qū)段的地表粗糙度后,運(yùn)用五種插值模型對(duì)各子區(qū)段的檢驗(yàn)樣本點(diǎn)進(jìn)行插值估算, 對(duì)比分析樣本點(diǎn)估算值與實(shí)測高程值, 得到測區(qū)不同地表粗糙度下插值模型的精度評(píng)價(jià)結(jié)果, 如圖3所示．

圖3 插值模型精度指標(biāo)對(duì)比Fig.3 Comparison of accuracy indicators of the interpolation methods in relation to surface roughness

由圖3可知, 各模型的插值誤差與河道地表粗糙度呈空間正相關(guān), 即地表粗糙度越高, 插值模型的精度越低,地表粗糙度對(duì)插值精度指標(biāo)的影響越明顯．當(dāng)Sr小于1.03時(shí), 各模型插值效果較為接近, 其中IDW的插值精度相對(duì)較高, TIN次之,此時(shí)地表粗糙度對(duì)插值模型精度的影響較小; 當(dāng)Sr為1.03～1.07時(shí), 各模型精度指標(biāo)的差距逐漸增大, IDW的精度優(yōu)勢不顯著,而TIN插值表現(xiàn)較為穩(wěn)健, MLS的精度變化較大; 當(dāng)Sr大于1.07時(shí), 各模型精度有所降低, 其中OK的精度最高, RBF次之, IDW和TIN最差．綜上得出, 當(dāng)Sr小于1.03時(shí), IDW的插值精度較高; 當(dāng)Sr為1.03～1.07時(shí), TIN的插值效果較好; 當(dāng)Sr大于1.07時(shí), OK的插值效果較好．

2.2 DEM插值結(jié)果

為了進(jìn)一步分析各插值模型在不同地形下的插值特點(diǎn), 分別選取地表粗糙度為1.011,1.048,1.083的三個(gè)典型子區(qū)段進(jìn)行DEM插值效果對(duì)比,插值結(jié)果如圖4所示．由圖4可得出,隨著地表粗糙度的增加,各模型對(duì)于地形走勢的預(yù)測也有所變化．Sr=1.011時(shí), 地表粗糙度較低, 地形起伏不明顯, 五種插值模型對(duì)于地形分布的預(yù)測結(jié)果較為接近;Sr=1.048時(shí), 五種模型的DEM插值結(jié)果差別顯著, 特別是局部區(qū)域的河谷變化,但河道輪廓的預(yù)測結(jié)果較為接近;Sr=1.083時(shí), 各模型插值效果受地形起伏變化影響更加顯著, 河道走勢預(yù)測差別較大, 其中OK和MLS的地形等高線分布較為規(guī)則, 而IDW, RBF和TIN的等高線則向岸邊擴(kuò)展, 近似呈月牙狀分布．

圖4 DEM插值結(jié)果Fig.4 Underwater DEM interpolation results

為了探究不同模型的插值特點(diǎn),結(jié)合各插值模型原理對(duì)上述研究結(jié)果進(jìn)一步分析．其中IDW作為一種精確的插值模型,插值時(shí)須保證樣本點(diǎn)分布的均勻性,而無人船在測量中無法始終保持勻速, 可能造成實(shí)測高程點(diǎn)分布不均的現(xiàn)象[10],進(jìn)而出現(xiàn)圖4(b)和圖4(c)中類似于“牛眼”的球狀突起(圖中圓圈位置),并且IDW只使用樣本點(diǎn)與估算點(diǎn)之間的距離作為權(quán)重指標(biāo),加權(quán)時(shí)對(duì)數(shù)據(jù)進(jìn)行平滑處理,會(huì)導(dǎo)致局部細(xì)節(jié)被忽略;TIN在任何地形環(huán)境下均使用平三角形來代替局部地形,對(duì)一些局部凹陷或淤泥抬升的區(qū)域會(huì)進(jìn)行平滑處理,進(jìn)而產(chǎn)生誤差,因此在地形變化平緩的區(qū)域插值精度較高,而在地形較為復(fù)雜的區(qū)域誤差較大;OK的空間適應(yīng)性強(qiáng)于其他插值模型,在選擇變異函數(shù)進(jìn)行插值時(shí),充分考慮測區(qū)整體地形的空間相關(guān)性,插值結(jié)果為最優(yōu)無偏估計(jì)量[11],能夠有效減小地形變化造成的估算誤差;RBF作為一種神經(jīng)網(wǎng)絡(luò)算法,雖然插值精度高、運(yùn)算速度快,但其對(duì)樣本點(diǎn)的依賴性較強(qiáng)[12],而無人船實(shí)測樣本點(diǎn)密度相對(duì)稀疏,無法對(duì)基函數(shù)進(jìn)行大量訓(xùn)練,進(jìn)而無法發(fā)揮算法優(yōu)勢;MLS是先擬合出光滑曲面,再對(duì)區(qū)域內(nèi)地形點(diǎn)進(jìn)行插值,擬合曲面由權(quán)函數(shù)控制,地形變化時(shí)不能自適應(yīng)地改變曲面形態(tài),預(yù)測精度受地形因素的影響最為顯著．

3 改進(jìn)算法

根據(jù)測區(qū)各插值模型的插值特點(diǎn),對(duì)插值流程進(jìn)行優(yōu)化設(shè)計(jì), 提出了一種能夠根據(jù)測量區(qū)域表面粗糙度自動(dòng)選擇插值方法的的改進(jìn)算法,主要步驟如下:

1) 將測區(qū)按照原始測線的布設(shè)狀況劃分為若干個(gè)區(qū)域面積接近,高程點(diǎn)分布均勻的子區(qū)段;

2) 依靠子區(qū)段原始測點(diǎn)的高程分布特征, 求取各子區(qū)段的地表粗糙度,根據(jù)地表粗糙度初步預(yù)測子區(qū)段水下地形的空間分布;

3) 判斷待插值點(diǎn)所在的子區(qū)段, 選取該子區(qū)段精度較高的插值模型對(duì)待插值點(diǎn)的高程進(jìn)行估算．

應(yīng)用此改進(jìn)算法對(duì)原始數(shù)據(jù)重新進(jìn)行插值計(jì)算, 得到三個(gè)典型子區(qū)段的DEM插值結(jié)果, 如圖5所示．

圖5 改進(jìn)算法DEM圖Fig.5 DEM of improved algorithm

為了驗(yàn)證改進(jìn)算法的精度, 在原始測區(qū)樣本集中重新抽取15%的高程點(diǎn)作為檢驗(yàn)樣本集,選擇各子區(qū)段的最優(yōu)插值模型,估算檢測樣本集,對(duì)比此算法與五種插值模型的預(yù)測精度, 結(jié)果如圖6所示．由圖6可知, 使用單一插值模型對(duì)測區(qū)整體檢驗(yàn)樣本估算時(shí),OK稍好,TIN次之,IDW和RBF誤差較大,MLS的預(yù)測吻合度最低; 改進(jìn)算法的η降低了15%左右,Re降低了10%左右,R2提高了0.005左右, 各項(xiàng)精度評(píng)價(jià)指標(biāo)均有所改善．從五種模型插值特點(diǎn)來看,使用某一種模型進(jìn)行樣本點(diǎn)高程估算時(shí),插值精度很大程度上依賴于局部地形特征,并且由于不同插值模型對(duì)地形變化敏感度的差異,選用單一模型會(huì)削弱樣本點(diǎn)的自然屬性,不能完整體現(xiàn)非線性特征數(shù)據(jù)的空間分布．而改進(jìn)算法能有效利用樣本點(diǎn)傳遞出的地形信息,充分發(fā)揮不同插值模型的優(yōu)勢,在一定程度上控制插值誤差,具有一定的實(shí)際工程適用性．

圖6 不同插值模型精度指標(biāo)對(duì)比Fig.6 Comparison of the accuracy indicators of different interpolation methods

4 結(jié)論

本文基于上海潘涇水下地形點(diǎn)實(shí)測數(shù)據(jù),選用五種插值模型,分析水下不同地形空間分布形態(tài)下的插值精度和插值效果,并提出一種改進(jìn)插值流程算法．得到如下結(jié)論:

1) 五種插值模型的測量誤差與地表粗糙度之間均呈空間正相關(guān),即隨著地形空間粗糙度的提高,各模型的插值精度均有下降;

2) IDW在測區(qū)地形變化不大的區(qū)域(Sr<1.03)插值精度相對(duì)較高, TIN在測區(qū)地表粗糙度為1.03～1.07時(shí)插值效果較為優(yōu)異, OK在測區(qū)地形變化復(fù)雜的區(qū)域(Sr>1.07)具有較高的插值精度和較強(qiáng)的穩(wěn)健性;

3) 改進(jìn)算法相較于傳統(tǒng)插值模型精度有所提高,具有較好的工程應(yīng)用前景．

本研究僅針對(duì)于潘涇地區(qū),在其他河流水系的可靠性有待驗(yàn)證．后續(xù)將在不同地域、不同水系開展無人船測量工作, 以探究插值模型在不同水下地形空間分布下的具體適用范圍．