蔡文婷
(鼎信信息科技有限責任公司,廣東廣州510623)
在創新科技時代引領下,以高新技術為支撐的海洋戰略將成為國家社會經濟發展的戰略要地[1]。“十二五”海洋科學和技術發展規劃指出[2]:“海洋調查需步入常態化和全球化,海洋觀測需進入立體觀測時代,并向實時化、系統化、信息化、數字化方向發展。”近海大陸架蘊含豐富的油氣、食物資源,優越的航運資源,是資源開發的戰略性基地。以資源開發為目的的各類工程建設,例如灘涂養殖、能源開發、海岸防護、港灣建設、航道開發、電纜鋪設等,都需要根據實際需求獲取合適比例尺的水下地形。
傳統船載聲學測量往往基于測深設備并結合定位設備進行探深測量。首先,分析測區特點進行航線規劃,生成合理的測深條帶網并進行測量,再利用數據插值手段進行數據加密,從而得到測區水下地形,其測量精度高但量測范圍有限[3-5]。隨著航空、航天遙感技術發展,通過多分辨率、多成像方式的遙感影像反演水深已成為水深測量領域的新途徑[6-9]。本文結合兩種探測結果,取長補短,構建測區大范圍水下地形數據。
本次數據獲取依托南方電網海南聯網路由檢測工程。海南聯網由中國南方電網責任有限公司投資建設,是國內第一個500 kV超高壓、長距離、大容量的跨海聯網工程。本次海底路由檢測測區位于瓊州海峽,海纜北側登陸點在廣東省湛江市徐聞縣南嶺村附近,南側登陸點在海南省海口市澄邁縣林詩島附近的玉包角。以基于多源數據大范圍無縫水下地形構建為目的,收集研究區多波束側掃聲納數據用于精確表征測區水下地形,收集全球30"海底地形數據與30 m全球DEM數據用于與側掃聲納數據融合構建適用于三維展示的大范圍水下地形。具體收集的資料如表1所示。

表1 研究區數據來源
考慮數據多源,為簡化運算、提高結果可用性,分別針對聲納數據及遙感數據先進行有針對性的預處理工作。
1.2.1 全球水下地形數據
數據空間分辨率為30",為提取水下地形數據,需遍歷數據柵格選取數值小于0的數據并將計算結果轉換為離散點。此外,還應剔除側掃聲納數據所覆蓋范圍。
1.2.2 全球Aster DEM數據
數據空間分辨率30 m,不包含水下地形數據,可直接將數據轉換為點數據。此外,由于數據空間分辨率高,可通過該數據提取海陸分界線,用于后續插值折斷線需要。
1.2.3 測區多波束側掃聲納數據
側掃聲納數據空間分辨率為1 m,直接參與插值則會因離散點密度較大造成與30"水下地形數據交界處的內插結果出現鋸齒狀。
為消除鋸齒現象,需對多波束側掃聲納數據進行離散點抽稀工作,具體算法如下:(1)構建分區統計格網;(2)遍歷格網內多波束點云,尋找最小值作為該格網中心水深值。
作為地質統計學的方法,克里金法已被廣泛應用于地質學、水文學、土壤學等研究“時空變量”的數據處理領域,特別是在測深數據處理方面[10]。
普通克里金插值顧及了各觀測點間的空間關系,可獲得未知變量的最優線性無偏估計[11],具有較高的估值計算精度。實際插值過程中,若存在異常數據,其變異函數會受異常值影響而偏離正確形狀,進而降低插值結果精度。
本文引入一種經驗貝葉斯克里金,實踐證明,該方法具有較好的抗差能力。
假設在待估計點(x)的海域附近共有n個實測點,即x1,x2,…,xn,其樣本水深值為Z(xi)。那么,普通克里金法的插值公式為:

式中,λi為權重系數,表示各空間樣本點xi處的水深觀測值Z(xi)對估計值Z*(x)的貢獻程度。
在引入變異函數的情況下,根據協方差與變異函數的關系γ(h)=c(0)-c(h),變異函數可等同于普通克里金方程組和克里金估計方差,即:

也可以將克立格方程組和估計方差用變異函數寫成上述矩陣形式:

在具體求解中,如何有效獲取變異函數的最優公式是克里金插值方法的關鍵,其主要受制于變異函數理論模型的選取和模型參數的估計。
2.1.1 樣本變異函數值計算
由于數據分布不規則,為使較多的數據點參加計算,使求解出的樣本變異函數值γ(dl)可信,應求解時取一常量Δd,尋找利用滿足的已知水深點對,距離dl的變異函數樣本值的求解公式為:

式中,ml為滿足的點對數目。
2.1.2 變異函數理論模型的擬合
理論變異函數模型較多,常用的有指數、高斯、線性、對數、塊金效應、冪函數、二次、推理二次、球狀和孔穴效應、立方和五球形[12]。
由于理論變異函數模型比較多,一般利用殘差、均方根預測誤差以及相對均方差等指標來衡量評定各種模型的有效性。
經驗貝葉斯克里金法與普通克里金方法有所不同,它通過估計基礎半變異函數來說明所引入的誤差。普通克里金方法由于不考慮半變異函數估計的不確定性,因此低估了預測的標準誤差。
經驗貝葉斯克里金由于可準確預測一般程度上不穩定的數據,對于小型數據集,比其他克里金法更準確。此外,該模型可以用于插入非平穩數據,以至于可以在較大區域內,局部地將數據改造為高斯分布。
經驗貝葉斯克里金法與普通克里金方法不同,它使用固有的0階隨機函數(IRF-0)作為克里金模型。
2.2.1 樣本變異函數模型
對于給定距離h,經驗貝葉斯克里金法使用以下形式的半變異函數:

塊金值和b(坡度)必須為正值,而α(冪)必須介于0.25和1.75之間。在EBK中,可以分析參數估計的經驗分布,因為在每個位置都估計了多個半變異函數。
2.2.2 半變異函數估計
與其他克里金法不同,EBK中的半變異函數參數由受限最大似然法估計。計算過程中,可以按以下方式估計半變異函數:(1)通過半方差圖模型對數據進行分析估計;(2)基于該半方差圖,生產每個輸入數據位置的新值;(3)使用新的模擬數據重新估計生成新的半方差圖,最后根據這個半方差圖的范圍計算需要使用的貝葉斯經驗規則。
為評價整體插值結果的優劣,本文根據以下四種方案評價本文所采用的技術路線的優勢:
方案一:不進行任何點云抽稀,以海陸分界線為折斷線,直接將所有數據參與運算;
方案二:對數據進行預處理工作,利用普通克里金方法進行插值;
方案三:選取已處理的數據,以海陸分界線為折斷線,利用普通克里金插值;
方案四:選取已處理的數據,以海陸分界線為折斷線,利用經驗貝葉斯克里金插值。
各方案插值結果如圖1所示。

圖1 各方案插值結果示意圖(相同配色方案)
從方案一插值結果可以看出,由于數據點云密度存在差異,多波束測量結果附近,由于搜索圓以最近距離為準則進行近鄰點搜索,造成近鄰點分布不均勻,插值結果存在明顯條帶階躍性;同理,由于后續方案采取點抽稀,近鄰點選取合理,插值結果過渡自然;從方案二配色方案可知,由于未利用海陸分界線作為折斷線,造成海陸分界處待插值點的近鄰點可能同時包含正負值,海陸分界錯誤。對比方案三、方案四,由于二者采取了點抽稀及折斷線等優化方式,整體插值效果較好,但分析交叉驗證結果可以發現,方案四的標準差較低,整體插值精度較高,如圖2所示。

圖2 兩種克里金插值交叉驗證誤差標準差
在現有測區多波束數據的基礎上結合全球30"海底地形數據與30 m全球DEM數據為數據源,通過數據整理與抽稀、海岸線折斷線以及經驗貝葉斯克里金插值,構建覆蓋測區的大范圍的無縫地形的技術流程可靠:在進行數據插值時,已知點的密度和搜索圓共同決定插值點相鄰點分布情況,在已知點的分布密度相差較大的情況下應進行相應的抽稀,以達到密度分布均勻的目的;涉及海陸分界的插值,為保證海岸線準確,應事先提取海岸線作為插值時的硬折斷線。本文所介紹的經驗貝葉斯克里金由于可準確預測一般程度上不穩定的數據,對于小型數據集,其插值結果比其他克里金法更準確。