基于多源數(shù)據(jù)的大范圍無(wú)縫水下地形構(gòu)建研究

蔡文婷

（鼎信信息科技有限責(zé)任公司，廣東廣州510623）

0 引言

在創(chuàng)新科技時(shí)代引領(lǐng)下，以高新技術(shù)為支撐的海洋戰(zhàn)略將成為國(guó)家社會(huì)經(jīng)濟(jì)發(fā)展的戰(zhàn)略要地[1]?！笆濉焙Ｑ罂茖W(xué)和技術(shù)發(fā)展規(guī)劃指出[2]：“海洋調(diào)查需步入常態(tài)化和全球化，海洋觀測(cè)需進(jìn)入立體觀測(cè)時(shí)代，并向?qū)崟r(shí)化、系統(tǒng)化、信息化、數(shù)字化方向發(fā)展?！苯４箨懠芴N(yùn)含豐富的油氣、食物資源，優(yōu)越的航運(yùn)資源，是資源開發(fā)的戰(zhàn)略性基地。以資源開發(fā)為目的的各類工程建設(shè)，例如灘涂養(yǎng)殖、能源開發(fā)、海岸防護(hù)、港灣建設(shè)、航道開發(fā)、電纜鋪設(shè)等，都需要根據(jù)實(shí)際需求獲取合適比例尺的水下地形。

傳統(tǒng)船載聲學(xué)測(cè)量往往基于測(cè)深設(shè)備并結(jié)合定位設(shè)備進(jìn)行探深測(cè)量。首先，分析測(cè)區(qū)特點(diǎn)進(jìn)行航線規(guī)劃，生成合理的測(cè)深條帶網(wǎng)并進(jìn)行測(cè)量，再利用數(shù)據(jù)插值手段進(jìn)行數(shù)據(jù)加密，從而得到測(cè)區(qū)水下地形，其測(cè)量精度高但量測(cè)范圍有限[3-5]。隨著航空、航天遙感技術(shù)發(fā)展，通過(guò)多分辨率、多成像方式的遙感影像反演水深已成為水深測(cè)量領(lǐng)域的新途徑[6-9]。本文結(jié)合兩種探測(cè)結(jié)果，取長(zhǎng)補(bǔ)短，構(gòu)建測(cè)區(qū)大范圍水下地形數(shù)據(jù)。

1 研究區(qū)與數(shù)據(jù)預(yù)處理

1.1 研究區(qū)簡(jiǎn)介

本次數(shù)據(jù)獲取依托南方電網(wǎng)海南聯(lián)網(wǎng)路由檢測(cè)工程。海南聯(lián)網(wǎng)由中國(guó)南方電網(wǎng)責(zé)任有限公司投資建設(shè)，是國(guó)內(nèi)第一個(gè)500 kV超高壓、長(zhǎng)距離、大容量的跨海聯(lián)網(wǎng)工程。本次海底路由檢測(cè)測(cè)區(qū)位于瓊州海峽，海纜北側(cè)登陸點(diǎn)在廣東省湛江市徐聞縣南嶺村附近，南側(cè)登陸點(diǎn)在海南省?？谑谐芜~縣林詩(shī)島附近的玉包角。以基于多源數(shù)據(jù)大范圍無(wú)縫水下地形構(gòu)建為目的，收集研究區(qū)多波束側(cè)掃聲納數(shù)據(jù)用于精確表征測(cè)區(qū)水下地形，收集全球30"海底地形數(shù)據(jù)與30 m全球DEM數(shù)據(jù)用于與側(cè)掃聲納數(shù)據(jù)融合構(gòu)建適用于三維展示的大范圍水下地形。具體收集的資料如表1所示。

表1 研究區(qū)數(shù)據(jù)來(lái)源

1.2 數(shù)據(jù)預(yù)處理

考慮數(shù)據(jù)多源，為簡(jiǎn)化運(yùn)算、提高結(jié)果可用性，分別針對(duì)聲納數(shù)據(jù)及遙感數(shù)據(jù)先進(jìn)行有針對(duì)性的預(yù)處理工作。

1.2.1 全球水下地形數(shù)據(jù)

數(shù)據(jù)空間分辨率為30"，為提取水下地形數(shù)據(jù)，需遍歷數(shù)據(jù)柵格選取數(shù)值小于0的數(shù)據(jù)并將計(jì)算結(jié)果轉(zhuǎn)換為離散點(diǎn)。此外，還應(yīng)剔除側(cè)掃聲納數(shù)據(jù)所覆蓋范圍。

1.2.2 全球Aster DEM數(shù)據(jù)

數(shù)據(jù)空間分辨率30 m，不包含水下地形數(shù)據(jù)，可直接將數(shù)據(jù)轉(zhuǎn)換為點(diǎn)數(shù)據(jù)。此外，由于數(shù)據(jù)空間分辨率高，可通過(guò)該數(shù)據(jù)提取海陸分界線，用于后續(xù)插值折斷線需要。

1.2.3 測(cè)區(qū)多波束側(cè)掃聲納數(shù)據(jù)

側(cè)掃聲納數(shù)據(jù)空間分辨率為1 m，直接參與插值則會(huì)因離散點(diǎn)密度較大造成與30"水下地形數(shù)據(jù)交界處的內(nèi)插結(jié)果出現(xiàn)鋸齒狀。

為消除鋸齒現(xiàn)象，需對(duì)多波束側(cè)掃聲納數(shù)據(jù)進(jìn)行離散點(diǎn)抽稀工作，具體算法如下：（1）構(gòu)建分區(qū)統(tǒng)計(jì)格網(wǎng)；（2）遍歷格網(wǎng)內(nèi)多波束點(diǎn)云，尋找最小值作為該格網(wǎng)中心水深值。

2 經(jīng)驗(yàn)貝葉斯克里金插值

作為地質(zhì)統(tǒng)計(jì)學(xué)的方法，克里金法已被廣泛應(yīng)用于地質(zhì)學(xué)、水文學(xué)、土壤學(xué)等研究“時(shí)空變量”的數(shù)據(jù)處理領(lǐng)域，特別是在測(cè)深數(shù)據(jù)處理方面[10]。

普通克里金插值顧及了各觀測(cè)點(diǎn)間的空間關(guān)系，可獲得未知變量的最優(yōu)線性無(wú)偏估計(jì)[11]，具有較高的估值計(jì)算精度。實(shí)際插值過(guò)程中，若存在異常數(shù)據(jù)，其變異函數(shù)會(huì)受異常值影響而偏離正確形狀，進(jìn)而降低插值結(jié)果精度。

本文引入一種經(jīng)驗(yàn)貝葉斯克里金，實(shí)踐證明，該方法具有較好的抗差能力。

2.1 普通克里金差值原理

假設(shè)在待估計(jì)點(diǎn)（x）的海域附近共有n個(gè)實(shí)測(cè)點(diǎn)，即x1，x2，…，xn，其樣本水深值為Z（xi）。那么，普通克里金法的插值公式為：

式中，λi為權(quán)重系數(shù)，表示各空間樣本點(diǎn)xi處的水深觀測(cè)值Z（xi）對(duì)估計(jì)值Z*（x）的貢獻(xiàn)程度。

在引入變異函數(shù)的情況下，根據(jù)協(xié)方差與變異函數(shù)的關(guān)系γ（h）=c（0）-c（h），變異函數(shù)可等同于普通克里金方程組和克里金估計(jì)方差，即：

也可以將克立格方程組和估計(jì)方差用變異函數(shù)寫成上述矩陣形式：

在具體求解中，如何有效獲取變異函數(shù)的最優(yōu)公式是克里金插值方法的關(guān)鍵，其主要受制于變異函數(shù)理論模型的選取和模型參數(shù)的估計(jì)。

2.1.1 樣本變異函數(shù)值計(jì)算

由于數(shù)據(jù)分布不規(guī)則，為使較多的數(shù)據(jù)點(diǎn)參加計(jì)算，使求解出的樣本變異函數(shù)值γ（dl）可信，應(yīng)求解時(shí)取一常量Δd，尋找利用滿足的已知水深點(diǎn)對(duì)，距離dl的變異函數(shù)樣本值的求解公式為：

式中，ml為滿足的點(diǎn)對(duì)數(shù)目。

2.1.2 變異函數(shù)理論模型的擬合

理論變異函數(shù)模型較多，常用的有指數(shù)、高斯、線性、對(duì)數(shù)、塊金效應(yīng)、冪函數(shù)、二次、推理二次、球狀和孔穴效應(yīng)、立方和五球形[12]。

由于理論變異函數(shù)模型比較多，一般利用殘差、均方根預(yù)測(cè)誤差以及相對(duì)均方差等指標(biāo)來(lái)衡量評(píng)定各種模型的有效性。

2.2 經(jīng)驗(yàn)貝葉斯克里金（EBK）

經(jīng)驗(yàn)貝葉斯克里金法與普通克里金方法有所不同，它通過(guò)估計(jì)基礎(chǔ)半變異函數(shù)來(lái)說(shuō)明所引入的誤差。普通克里金方法由于不考慮半變異函數(shù)估計(jì)的不確定性，因此低估了預(yù)測(cè)的標(biāo)準(zhǔn)誤差。

經(jīng)驗(yàn)貝葉斯克里金由于可準(zhǔn)確預(yù)測(cè)一般程度上不穩(wěn)定的數(shù)據(jù)，對(duì)于小型數(shù)據(jù)集，比其他克里金法更準(zhǔn)確。此外，該模型可以用于插入非平穩(wěn)數(shù)據(jù)，以至于可以在較大區(qū)域內(nèi)，局部地將數(shù)據(jù)改造為高斯分布。

經(jīng)驗(yàn)貝葉斯克里金法與普通克里金方法不同，它使用固有的0階隨機(jī)函數(shù)（IRF-0）作為克里金模型。

2.2.1 樣本變異函數(shù)模型

對(duì)于給定距離h，經(jīng)驗(yàn)貝葉斯克里金法使用以下形式的半變異函數(shù)：

塊金值和b（坡度）必須為正值，而α（冪）必須介于0.25和1.75之間。在EBK中，可以分析參數(shù)估計(jì)的經(jīng)驗(yàn)分布，因?yàn)樵诿總€(gè)位置都估計(jì)了多個(gè)半變異函數(shù)。

2.2.2 半變異函數(shù)估計(jì)

與其他克里金法不同，EBK中的半變異函數(shù)參數(shù)由受限最大似然法估計(jì)。計(jì)算過(guò)程中，可以按以下方式估計(jì)半變異函數(shù)：（1）通過(guò)半方差圖模型對(duì)數(shù)據(jù)進(jìn)行分析估計(jì)；（2）基于該半方差圖，生產(chǎn)每個(gè)輸入數(shù)據(jù)位置的新值；（3）使用新的模擬數(shù)據(jù)重新估計(jì)生成新的半方差圖，最后根據(jù)這個(gè)半方差圖的范圍計(jì)算需要使用的貝葉斯經(jīng)驗(yàn)規(guī)則。

3 實(shí)驗(yàn)結(jié)果及精度評(píng)價(jià)

為評(píng)價(jià)整體插值結(jié)果的優(yōu)劣，本文根據(jù)以下四種方案評(píng)價(jià)本文所采用的技術(shù)路線的優(yōu)勢(shì)：

方案一：不進(jìn)行任何點(diǎn)云抽稀，以海陸分界線為折斷線，直接將所有數(shù)據(jù)參與運(yùn)算；

方案二：對(duì)數(shù)據(jù)進(jìn)行預(yù)處理工作，利用普通克里金方法進(jìn)行插值；

方案三：選取已處理的數(shù)據(jù)，以海陸分界線為折斷線，利用普通克里金插值；

方案四：選取已處理的數(shù)據(jù)，以海陸分界線為折斷線，利用經(jīng)驗(yàn)貝葉斯克里金插值。

各方案插值結(jié)果如圖1所示。

圖1 各方案插值結(jié)果示意圖（相同配色方案）

從方案一插值結(jié)果可以看出，由于數(shù)據(jù)點(diǎn)云密度存在差異，多波束測(cè)量結(jié)果附近，由于搜索圓以最近距離為準(zhǔn)則進(jìn)行近鄰點(diǎn)搜索，造成近鄰點(diǎn)分布不均勻，插值結(jié)果存在明顯條帶階躍性；同理，由于后續(xù)方案采取點(diǎn)抽稀，近鄰點(diǎn)選取合理，插值結(jié)果過(guò)渡自然；從方案二配色方案可知，由于未利用海陸分界線作為折斷線，造成海陸分界處待插值點(diǎn)的近鄰點(diǎn)可能同時(shí)包含正負(fù)值，海陸分界錯(cuò)誤。對(duì)比方案三、方案四，由于二者采取了點(diǎn)抽稀及折斷線等優(yōu)化方式，整體插值效果較好，但分析交叉驗(yàn)證結(jié)果可以發(fā)現(xiàn)，方案四的標(biāo)準(zhǔn)差較低，整體插值精度較高，如圖2所示。

圖2 兩種克里金插值交叉驗(yàn)證誤差標(biāo)準(zhǔn)差

4 結(jié)語(yǔ)

在現(xiàn)有測(cè)區(qū)多波束數(shù)據(jù)的基礎(chǔ)上結(jié)合全球30"海底地形數(shù)據(jù)與30 m全球DEM數(shù)據(jù)為數(shù)據(jù)源，通過(guò)數(shù)據(jù)整理與抽稀、海岸線折斷線以及經(jīng)驗(yàn)貝葉斯克里金插值，構(gòu)建覆蓋測(cè)區(qū)的大范圍的無(wú)縫地形的技術(shù)流程可靠：在進(jìn)行數(shù)據(jù)插值時(shí)，已知點(diǎn)的密度和搜索圓共同決定插值點(diǎn)相鄰點(diǎn)分布情況，在已知點(diǎn)的分布密度相差較大的情況下應(yīng)進(jìn)行相應(yīng)的抽稀，以達(dá)到密度分布均勻的目的；涉及海陸分界的插值，為保證海岸線準(zhǔn)確，應(yīng)事先提取海岸線作為插值時(shí)的硬折斷線。本文所介紹的經(jīng)驗(yàn)貝葉斯克里金由于可準(zhǔn)確預(yù)測(cè)一般程度上不穩(wěn)定的數(shù)據(jù)，對(duì)于小型數(shù)據(jù)集，其插值結(jié)果比其他克里金法更準(zhǔn)確。