MF多源測深數據融合方法及大洋水深模型構建

劉 洋，吳自銀，趙荻能，周潔瓊，尚繼宏，王明偉,朱 超，魯號號

1. 自然資源部第二海洋研究所，浙江 杭州 310012； 2. 國家海洋局海底科學重點實驗室，浙江 杭州 310012； 3. 上海交通大學海洋研究院，上海 200030

精確可靠的水深數據幾乎是從事一切海洋活動的基礎[1]，高分辨率的數字水深模型(digital bathymetric model，DBM)不僅可用于海圖的構建，還可用于地貌特征與成因[2-3]、海底構造、海洋資源、水動力流體模型、生物棲息地與生態系統等多個方面的研究[4-14]。對海洋測深的系統研究持續了一個多世紀[15-19]，從早期極低效率的錘線測深發展到當前高精度、高分辨率、高效率的多波束測深技術[20-22]。然而，由于全球海洋極其廣闊，基于船載模式的多波束測深成本高昂[23]，導致目前仍有超過80%的全球海域未被多波束探測[24]。因此，在構建深海大范圍DBM時仍需使用歷史測深資料[25]。當前大洋測深數據主要針對不同目的獨立采集，導致來源、格式和分辨率等方面存在諸多差異[26-28]。全球大洋測深數據的多源復雜性，不利于高精度DBM的構建，更不利于日益迫切認識海洋、經略海洋的國家戰略需求。

DBM的構建也稱數據的網格化。網格化是二維采樣、插值和二維外推的結合[29]。如何基于插值算法將多源水深數據融合進而構建DBM一直是學者們的研究熱點。根據目前國內外對多源水深數據融合技術的研究，融合方法主要是基于克里金、雙三次樣條、張力樣條等插值算法進行綜合應用，如文獻[30]將多波束、單波束、歷史海圖等水深數據進行融合，編制了北冰洋的國際海底地形圖(international bathymetric chart of the arctic ocean，IBCAO)；文獻[31]根據多種水深數據匯編了南大洋的國際海底地形圖(international bathymetric chart of the southern ocean，IBCSO)；文獻[32]結合不同時期的多波束、單波束和海岸線等數據，構建了南極洲東部George V大陸邊緣的DBM，得到該區較詳盡的海底地貌信息。上述方法雖實現了多源水深數據的融合，但融合后的數據空間分辨率低(>500 m)，且難以保證數據稀疏區和空白區的插值精度。因此，在保證數據稀疏區插值精度的同時，既要提高數據空白區的精度，又要保留高分辨率水深源數據的細節特征，是值得深入研究的科學技術問題。

為解決上述問題，本文提出合并-融合(merge-fusion，MF)法。該方法將多波束、單波束、電子海圖數據與通用全球海洋地形(general bathymetric chart of the oceans，GEBCO)數據融合，在保留高分辨率地形細節特征的同時，合理填補了數據空白區，形成“挑戰者深淵”分辨率100 m的高精度DBM。結果表明，該方法融合的DBM能更好地反映精細的地形特征信息，為區域性或全球性DBM的構建以及多源異構水深數據的管理利用提供了實用化技術方案。

1 MF方法基本原理與步驟

MF(merge-fusion)法的核心思想為“合并-融合”。如圖1所示，該方法主要包括6個步驟：

(1) 原始數據處理與入庫。根據不同水深數據的來源、類型，對其進行精細處理。以多波束數據為例，對于系統誤差，通過潮位改正、聲速剖面改正、換能器安裝偏差校準等方法進行消除[33-39]；對于異常值，通過人工或者自動濾波算法識別并剔除[20]；而隨機誤差，只能通過對測深數據的精度評估評價其誤差大小[40]。建立多源水深數據庫，對清理后的水深數據進行統一管理。

(2) 多源數據精度評估。對多源水深數據進行精度評估是保證成果網格數據質量的重要手段[40]。在同一水域，以多波束條帶數據的中央波束測深結果為基準，使用平均差值式[30](1)對多源數據進行精度評估

(1)

式中，δ為平均水深差；q為測深數據個數；DMi為多波束水深數據；DSi為多源水深數據，水深值取負值。

根據《GB 12327—1998 海道測量規范》，將水深偏差與水深的2%比較，若偏差小于水深的2%，符合規范要求不做處理；否則根據偏差值以多波束條帶數據的中央波束測深結果為基準，對其他類型數據進行最小二乘法平差處理[28,41]。通過式(2)來表示誤差效應

F(x,y)=a0+a1x+a2y+a3x2+a4xy+a5y2

(2)

式中，(x,y)為測點坐標；a0,a1,a2,a3,a4,a5為待定系數。

設水深觀測值Z的改正數為V，以多波束條帶數據的中央波束測深結果為基準，建立如下誤差方程

-(VSi-VMi)=(a0+a1x+a2y+a3x2+

a4xy+a5y2)Si-(a0+a1x+

a2y+a3x2+a4xy+a5y2)Mi-

(ZSi-ZMi)

(3)

式中，VSi、VMi分別為多源數據和多波束數據在點i位置(x,y)處的改正數；ZSi，ZMi分別為相應位置處的多源數據和多波束數據水深值。

圖1 MF方法原理及流程Fig.1 Principle and workflow of the MF method

根據∑(VSi-VMi)2=min解算誤差方程(3)得到誤差改正系數，最后依據式(2)對所有測深數據進行改正，即獲得平差后水深值。

(3) 多源數據網格化。為避免網格內潛在異常值的影響[42]，對多源水深數據庫導出的離散水深點進行中值濾波。中值濾波首先確定奇數像元網格G，網格內各像元按數值大小排序，并將網格節點內的中值水深d(x,y)作為原數值s(x,y)[43]

d(x,y)=median{s(x-u,y-v),(u,v)∈G}

(4)

為保證數據稀疏區域插值的準確性[44]，對濾波后的多源水深數據進行網格化[45-47]，形成低分辨率網格。地球科學領域最常見的網格化方法之一是雙三次樣條插值算法[48]，該算法可通過已知點精確插值，得到光滑的插值曲面，但在地形突變處或數據稀疏區易產生異常值，掩蓋真實地形。本文采用目前國際上普遍使用的連續曲率張力樣條插值算法[48-49]對多源數據進行網格化，該算法由于引入張力因子，可消除最小曲率擬合面中存在的較大波動和無關變形點[43]，網格化后的結果更加平滑、準確。該算法基于最小曲率擬合面原則

(5)

式中，(xi,yi,zi)為已知點坐標；N為已知點個數；fi為點(xi,yi)處格林函數線性組合解中的系數；δ(x-xi,y-yi)為給定響應函數；2為拉普拉斯算子。

邊界條件為

(6)

在此基礎上，引入張力因子T，得到

(1-TI)2(2z)-TI2z=

(7)

式中，TI為擬合面內部張力因子，在[0,1)區間內取值。

邊界條件為

(8)

式中，TB為擬合面邊界處張力因子，在[0,1)區間內取值。

當張力因子T=0時，式(7)簡化為式(5)，此時無張力，得到最小曲率平滑解；T=1時，張力最大，得到調和樣條解。對于大洋水深數據，張力因子一般選取0.32～0.4之間[30]。隨后對低分辨率網格濾波，將濾波后的低分辨率網格重采樣至目標分辨率，得到“基礎水深網格”。

(4) 多分辨率網格合并。為合并多分辨率網格，且保留高精度測深區域地形細節特征，使用高分辨率網格水深值替換低分辨率網格水深值。首先，將預處理步驟中輸出的高分辨率網格與基礎水深網格疊加分析，得到二者在同一位置處z的差值；對z的差值進行判斷，若差值非零則用多波束數據網格點內相應值進行替換，得到“合并水深網格”。

為消除網格合并過程出現的邊緣效應[50]，在高分辨率網格數據一側定義過渡帶，過渡帶寬度至少為3列網格點，具體范圍根據接邊情況來確定。其原理基于雙曲線加權函數

(9)

式中，L代表到下一高分辨率數據約束單元或過渡帶外約束單元的距離；w表示權重[31]。

由式(9)可看出，距離高分辨率數據網格單元越近，受深度影響越小，由此弱化邊緣效應發生的概率。

(5) 數據空白區融合。在數據空白區域范圍內，將GEBCO網格與合并水深網格“作差”，即GEBCO網格內每個數據點減去相同位置處插值數據點的值，得到包含經度、緯度和水深差異的差值文件。隨后網格化中值濾波后的差值文件，得到差值網格。將差值網格疊加整合到合并水深網格，得到融合后的DBM。

(6) 網格可視化及清理。經上述步驟得到的DBM難免存在各種“毛刺”等異常點，使用三維地形瀏覽法找出這些問題。為突出數據噪聲，可使用光照效果加以突出，最后通過曲面擬合算法[51]對網格問題區域進行編輯和誤差處理，以獲得最佳的數字水深模型。曲面擬合算法的構建原理：根據波束點擬合海底趨勢面，計算實測水深數據與趨勢面間的深度差值，結合誤差處理理論剔除異常點[52]。曲面擬合法中擬合函數的一般形式

(x,y,z)∈Q(xc,yc,zc)}

(10)

式中，Ψ(x,y)為曲面擬合函數；(x,y,z)為波束點空間坐標；hlm為多項式系數；k為多項式總階數；(xc,yc,zc)為待檢測點c；Q(xc,yc,zc)為以c為中心的曲面擬合函數的局部曲面擬合范圍。

2 數據來源

位于馬里亞納海溝南端的“挑戰者深淵”是全球海洋最深處，被譽為“世界第四極”，一直以來都是全球海洋科學研究的熱點區域[53-56]。近幾年來，我國在該區域組織了數次綜合海洋調查。高精度的海溝區DBM是開展載人深潛、海底原位觀測和板塊俯沖數值模擬[57]的基礎，但該區現有數據分辨率低、精度差，難以滿足海洋科學研究和工程實施要求。同時，該海區數據來源多樣，且存在諸多數據空白區，適合數據融合技術的方法研究。

本文數據主要來自于公開的海洋數據，包括多波束數據、單波束數據、電子海圖數據以及GEBCO數據。

2.1 多波束數據

構建此次DBM的核心數據主要是多波束數據，采集自2003年以后，且精度較高。在該區域共收集25個航次的多波束數據，其中，15個調查航次數據源于美國大氣與海洋局(NOAA)下屬的國家地質與地球物理數據中心(NGDC)，另外10個多波束測深航次數據則從日本海洋-地球科學和技術研究中心(JAMSTEC)獲得。

2.2 單波束數據

除多波束數據外，還從NGDC收集了2000年以來的大量單波束測深點數據(圖2)。由于收集的單波束數據未提供測深儀及定位系統等相關詳細信息，很難對所有數據進行精度評估。且單波束測深點空間分布非常不均勻，只在使用多波束探測區域分布較為密集，例如美國關島附近。

2.3 電子海圖數據

從美國國家海洋和大氣管理局(NOAA)下的海岸調查部門(Office of Coast Survey)共收集該區域2005年來的4幅電子海圖數據。除電子海圖數據外無相關元數據信息，因此無法全面評估其數據精度。

G區為多源數據精度評定所選區圖2 “挑戰者深淵”位置和數據源Fig.2 Location and data sources of the “Challenger Deep”

2.4 GEBCO數據

雖然已有大量水深數據，但仍存在范圍達上千平方公里的數據空白區，本文使用GEBCO_2014[58]水深數據集對數據空白區進行填補。GEBCO數據是當前國際主流水深數據集之一，而GEBCO_2014為目前全球測深數據集的最新版本，主要源于實測數據和利用衛星高度計反演的預測水深，網格分辨率為30″。

3 MF方法的實現

根據上述方法對實驗區多源水深數據進行融合處理，構建馬里亞納海溝“挑戰者深淵”的DBM。首先將多源水深數據處理為離散數據格式，轉換水平坐標，投影方式和深度基準面分別統一至UTM投影和平均海平面。隨后進行系統誤差處理，并使用CUBE算法清理多源水深數據的異常點[59-60]，后重新導出為離散格式，形成“挑戰者深淵”的多源水深數據庫。

為對多源水深數據進行精度評定，比較了美國關島西北部區域(圖2區域G)內多源水深數據的深度差異。該區域被多波束數據全覆蓋，精度較高，適合用于精度評估。離散的電子海圖和單波束水深點分別與分辨率100 m的多波束條帶數據的中央波束測深結果進行比較(圖3)。兩者分別與多波束進行一一組合，組合的準則是二者間距小于50 m[30]。選取115對電子海圖水深對比數據、207對單波束水深對比數據，依據式(1)得到對比結果，分別見表1和圖3。

藍色點代表電子海圖數據，紅色點代表單波束數據，直方圖表示深度差值絕對值分布圖3 電子海圖、單波束與多波束數據水深值對比Fig.3 Comparison of water depth of electronic charts, single-beam and multi-beam bathymetric data

Tab.1 Comparison of multi-source dataaccuracyassessmentm

對象平均差值平均水深±2%·水深電子海圖vs多波束-18.9-1039±20.7單波束vs多波束-23.2-1278±25.5GEBCO_2014vs多波束-10.7-1185±23.7

根據《GB 12327—1998 海道測量規范》，對于電子海圖數據，平均水深1039 m時，平均水深差值優于水深值的2%；對于單波束數據，平均水深1278 m時，平均水深差值同樣優于水深的2%。經過評定，單波束數據、電子海圖數據精度在限差容許范圍內，可作為本次數據融合實驗的數據集。

以區域G內多波束條帶數據的中央波束測深結果為基準，對GEBCO_2014數據進行精度評估。在分辨率100 m基礎上，得到二者水深差值的空間位置分布(圖4)。應用式(1)從29605對GEBCO對比數據中，得到對比結果，分別見表1和餅狀圖(圖4)。平均水深1185 m時，平均深度差優于水深的2%。GEBCO_2014數據精度同樣在限差容許范圍內，可作為本次數據空白區融合實驗的數據集。

圖4 GEBCO_2014數據與多波束數據水深值對比Fig.4 Comparison of water depth of GEBCO_2014 data and multi-beam bathymetric data

對多源水深數據(圖5(a))進行網格大小1000 m的中值濾波。隨后網格化濾波后的水深數據，形成分辨率1000 m的低分辨率水深網格(圖5(b))。對低分辨率網格進行平滑和濾波，將濾波后的低分辨率網格重采樣至100 m作為“基礎水深網格”(圖5(c))。

通過對多波束數據進行處理，例如聲速剖面改正、異常值剔除等測深條帶拼接步驟[33,43,61-63]得到分辨率100 m的高分辨水深網格(圖5(d))。然后將高分辨率水深網格與基礎水深網格疊加分析，得到二者在同一位置處z的差值。若某點處z的差值非零，則用高分辨率水深網格內相應水深值替換，得到“合并水深網格”(圖5(e))。為消除多分辨率網格間的邊緣效應，在靠近高分辨率網格數據邊界定義5 km的過渡帶(圖6)，極大地降低了邊緣效應產生的概率。

圖5 多源數據融合過程示例Fig.5 An example of multi-source data fusion

圖6 過渡帶定義前后對比Fig.6 Comparison of the DBM

結合GEBCO_2014數據，對數據空白區域進行融合。首先將GEBCO_2014網格與合并水深格網疊加合并，得到數據空白范圍內同一位置處的水深差值文件，差值文件包括經度、緯度和深度差值信息。隨后對差值文件進行網格大小1000 m的中值濾波。將濾波后的差值文件網格化為分辨率100 m的差值網格。最后，將差值網格與合并水深網格疊加整合，得到融合后的DBM(圖5(f))。

可視化瀏覽合并-融合后的水深網格，識別并高亮異常點，對網格問題區域進行編輯和誤差處理，最終得到滿意的DBM。

4 分析與討論

4.1 數據拼接對比

圖7為使用MF法得到的馬里亞納海溝“挑戰者深淵”DBM，網格分辨率100 m。經統計，本次所選區域多波束、單波束和電子海圖數據的覆蓋率分別為42.2%、0.5%、0.2%，使用GEBCO數據對剩余57%的數據空白區進行填補。通過網格可看出，MF方法不僅保留了高分辨率數據區詳細的海底特征，還為數據稀疏區和空白區提供了平滑的地形表面。

在馬里亞納海溝“挑戰者深淵”地區，除本文構建的“挑戰者深淵”DBM之外，只有全球數據集涵蓋整個大洋。本文將“挑戰者深淵”DBM與GEBCO_2014DBM進行比較。結果表明，在數據接邊處，MF法取得較為理想的效果，接邊處總體限差符合《GB 12327—1998 海道測量規范》，即水深大于100 m時，深度測量的極限誤差應優于水深的2%。在實驗區選取接邊差異較大區域M(圖7)，針對MF法和傳統方法得到剖面a-a′(圖8(a))和剖面b-b′(圖8(b))。從剖面圖可看出，兩種方法形成差異主要集中在距剖面線起點2～4 km范圍內區域。以剖面線起點為原點，剖面方向為正方向，選取(2.6 km，2.8 km)范圍內剖面A和B(圖8)，在水深約3600 m處，水深限差50 m時，對于傳統方法，平均坡度為72.32°(圖8(b))；而MF法，平均坡度為37.63°，與前者相比整體過渡較為平緩(圖8(a))。

4.2 數據細節對比

高分辨率水深數據結合越多，網格分辨率就會越高，這使得本次融合的分辨率100 m的“挑戰者深淵”DBM成為該區域目前精度最高的數據集；同時，該數據集對“挑戰者深淵”地區的地形細節特征刻畫有著顯著改善。例如，在網格單元100 m的“挑戰者深淵”DBM中，可清晰地看到海底峽谷等地形特征(圖9(d))，以及形態、輪廓更明顯的海山(圖9(c))，而在GEBCO_2014數據集中幾乎看不到此類地形特征。在數據稀疏區，利用清理過的單波束或電子海圖數據形成的低分辨率網格避免了異常點的出現(圖9(a)和9(b))。但仍存在數據密度和質量較差地區，這些區域包含潛在異常點、噪聲點，需進一步調查才能識別和處理。

5 結 論

(1) 提出并實現了一種適用于多源水深數據融合的MF方法，并給出詳細的融合技術方法流程。其基本原理是基于多源水深數據形成兩種分辨率網格，然后將多分辨率網格合并，在此基礎上使用GEBCO數據進行融合，疊加、合并得到最終的DBM。

M為數據拼接對比區；a、b、c、d 為數據細節對比區圖7 “挑戰者深淵”DBMFig.7 The DBM of the “Challenger Deep”

圖8 數據接邊處比較Fig.8 Comparison of data transition area

a1、b1、c1、d1為MF法所得DBM，a2、b2、c2、d2為傳統方法所得DBM，A代表異常點，C代表海底峽谷圖9 “挑戰者深淵”DBM與GEBCO_2014DBM細節對比Fig.9 Comparison of data details between “Challenger Deep” DBM and GEBCO_2014 DBM

(2) 基于馬里亞納海溝“挑戰者深淵”地區公開的多波束、單波束、電子海圖等水深數據，采用MF多源數據融合方法，構建了“挑戰者深淵”的高精度DBM，驗證了本文方法的可靠性和實用價值。

(3) 針對同一區域，利用本文方法得到的DBM與目前國際主流的GEBCO_2014數據集進行對比，表明本文方法更能夠體現高分辨率水深數據和數據空白區的細節信息、防止邊緣效應和異常點的產生，從而獲得高精度的海底DBM。