本文内容來源于《測繪學報》2024年第2期（審圖号GS京(2024)0297号）
利用擾動重力資料反演海底地形的高斯曲面函數估計方法
翟振和
, 孫中苗, 管斌, 馬健, 李端
地理資訊工程國家重點實驗室, 北京 100029
基金項目：國家自然科學基金（42174001；41774018）；地理資訊工程國家重點實驗室基金（SKLGIE2023-22-5）
摘要：利用2017年1月—2020年12月的SARAL衛星測高資料反演得到南海海域1'×1'格網擾動重力資料，通過與船載重力資料比較，精度達到5.5 mGal。提出了利用重力資料及高斯曲面函數求解區域特征參數進而估計海底地形模型的方法，利用SARAL衛星反演獲得的重力資料在南海海域開展了計算試驗。結果表明，利用ETOPO-1先驗模型，在10×10格網一組劃分條件下估計得到的1'×1'海底地形精度相對于先驗模型提高了約10 m；利用DTU18先驗模型，在9×9格網一組劃分條件下估計得到的1'×1'海底地形精度相對于先驗模型提高了約9 m。計算結果表明，利用衛星測高反演得到的重力資料，通過高斯曲面函數解算獲得的5個特征參數，可以在一定程度上代表相應區域海底地形的曲面特征，進而可在不依賴外部實測資料條件下對先驗海底模型進行精化。對于曲面估計而言，格網劃分越小，曲面函數就越能反映區域變化特征。是以，對于未來衛星測高技術發展而言，更高分辨率的重力場探測技術有望繼續提升海底地形細節的反演能力。
關鍵詞：擾動重力 海底地形 SARAL衛星 高斯曲面估計 ETOPO-1 DTU18
翟振和, 孫中苗, 管斌, 等. 利用擾動重力資料反演海底地形的高斯曲面函數估計方法[J]. 測繪學報，2024，53(2)：231-238. DOI: 10.11947/j.AGCS.2024.20230204ZHAI Zhenhe, SUN Zhongmiao, GUAN Bin, et al. An estimation method of seabed topography based on Gauss surface function using ocean gravity data[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(2): 231-238. DOI: 10.11947/j.AGCS.2024.20230204 閱讀全文：http://xb.chinasmp.com/article/2024/1001-1595/20240203.htm
引 言
海底地形資訊對于船艇航行、地質資源勘查及地球科學研究等具有重要意義，目前，海底地形資訊的擷取主要采用船載多波束測量、機載雷射測深、多光譜反演及衛星測高反演等技術[1-6]。船載多波束測量精度最高，已經成為海域高精度海底地形測繪的主要手段，但受限于測量平台和海洋客觀環境條件，船載測量技術很難達到全球均勻覆寫。海洋衛星測高通過雙星串飛、合成孔徑高度計等手段可在2.5年左右時間擷取全球高達1′分辨率的重力資料[7-9]，進而可反演獲得一定分辨率且均勻分布的海底地形，現已成為海底地形三維模型建構的重要手段之一。大陸自主研制的新一代海洋測高衛星于2023年3月發射，有望進一步提升寬闊海域海底地形資訊的精細程度。目前，各國研究機構廣泛應用包括船載測深、衛星測高等資料建構海底地形模型，美國國家環境資訊中心(National Centers for Environmental Information，NCEI)、丹麥科技大學(Technical University of Denmark，DTU)及斯克裡普斯海洋研究所(Scripps Institution of Oceanography，SIO)等多個研究機構相繼釋出了ETOPO-1、DTU18、SIO V19.1等全球海底地形格網模型[10-11]。國内學者基于衛星測高資料和船載水深資料也獲得了局部或全球海域海底地形，如BAT_WHU2020[12]。在海底地形反演方法上，目前廣泛利用衛星測高獲得的重力異常資料通過重力地質法、解析法、頻域法及最小二乘配置法等反演海底地形[13-25]，上述方法一般采用線性化後取1階項進行計算，文獻[26]顧及高階項影響建構了非線性2次相關模型，利用該模型可在已知少量船載水深/重力資料條件下反演得到一定精度的海底地形。盡管利用衛星測高資料反演海底地形取得了很大進展，但考慮到重力場與海底地形之間的複雜函數關系，仍有必要加大研究并拓展反演思路。本文基于重力場與海底地形的随機相關特點，提出通過建構曲面特征的思路在先驗模型基礎上進行推估得到新的海底地形，為進一步拓展衛星測高資料的應用價值提供借鑒。
1 考慮奇異影響的海洋擾動重力反演模型
利用衛星測高資料反演重力場的一般方法是在垂線偏差基礎上反演得到重力異常，考慮到擾動重力在大地邊值問題上的潛在優勢，本文采用擾動重力精确快速反演方法[27]計算得到格網點的擾動重力δgp (1)
式中，γ表示正常重力；ξ、η分别表示垂線偏差的子午分量和卯酉分量；α表示計算點到觀測點的方位角，核函數K(ψpq)的導數表示為
(2)
當球面角距ψ趨近0時，式(2)會出現奇異問題，為了避免奇異，仿照文獻[2]的思路在忽略二階及以上項影響條件下推導得到擾動重力在奇異點處的計算公式為
(3)
式中，ξy、ηx分别表示垂線偏差子午分量在緯度方向的梯度、垂線偏差卯酉分量在經度方向的梯度；s0表示格網半徑的大小。
2 海域重力資料反演海底地形的高斯曲面函數估計方法
在已有學者研究中可以發現，在局部區域重力擾動/異常與地形起伏之間有較強相關性，為此本文提出這樣一種假設，即在小區域重力資料的曲面特征與對應區域海底地形的曲面特征具有較強一緻性，如果已知精确的重力場特征資訊，就有可能恢複出相應的海底地形特征，并以此可估計出小區域的海底地形。為了對曲面變化特征進行定量分析，本文采用高斯曲面函數進行描述。二維高斯曲面函數可以用來描述兩個次元的資料分布情況，在圖像濾波、邊緣檢測等方面發揮着重要作用，可以生成高斯卷積核，用于圖像進行中的高斯濾波，可有效去除高斯噪聲，在機器學習核資料挖掘中也能發揮分類、聚類和回歸等作用，進而更加準确地預測資料。高斯曲面函數具體表現為5個特征參數即區域峰值、區域峰值對應的經度、緯度，以及區域資料在經度和緯度方向上的變化幅度。
假設衛星測高反演或其他手段獲得的區域格網重力資料為δg(i)(本文設定格網大小為1′×1′)，i=1，2, …, n。将資料按照從北至南、從西至東劃分為由m×m個1′×1′格網組成的小區域(m>3)。在每個小區域，利用已知重力資料按照高斯曲面函數模型(式(4))及最小二乘原理求解得到5個特征參數GT、GxT、GyT、CGx、CGy
(4)
式中，GT表示局部區域内重力場(如重力異常、擾動重力)的峰值；GxT、GyT表示峰值對應的經度、緯度資訊；CGx、CGy表示經度、緯度方向重力場變化的幅度；xi、yi表示每個重力場變量對應的經緯度資訊。
引入先驗海底地形模型Hs(i)，i=1，2, …, n，格網大小為1′×1′。利用先驗海底地形資料按照高斯曲面函數模型(式(5))及最小二乘原理求解得到5個特征參數HT、HxT、HyT、CHx、CHy
(5)
式中，HT表示局部區域内海底地形的峰值；HxT、HyT表示峰值對應的經度、緯度資訊；CHx、CHy表示經度、緯度方向海底地形變化的幅度；xi、yi表示每個海底地形變量對應的經緯度資訊。
從已有陸地重力、高程及海洋重力、水深資料(含格網資料和測線資料)分析來看，兩種資料在變化的峰值位置是基本一緻的，是以可将GxT、GyT表示待估計海底地形峰值對應的經度、緯度，将CGx、CGy賦以一定比例因子(式(6)、式(7))，而後作為待估計海底地形在小區域的變化幅度
(7)
将待估計海底地形的峰值HTnew在(HT-1000，HT+1000)範圍内循環，将scale在(-1/10，1/10)範圍内循環，上述循環為雙重循環同時進行，循環次序對結果沒有影響。而後按照式(8)進行求解
(8)
當小區域估計海底地形Hest－i和先驗海底地形資料Hs(i)之間的均方差小于一定門檻值時，循環結束，輸出小區域估計海底地形。
重複以上步驟，可将所有小區域計算完畢，最終輸出得到完整區域的海底地形，整個估計方法的計算流程如圖 1所示。
| |
| 圖 1 海底地形反演計算流程Fig. 1 Flowchart of seabed topography inversion |
| 圖選項 |
3 計算試驗
3.1 海域擾動重力反演
試驗資料采用法國空間局釋出的SARAL衛星GDR-F版資料，時間跨度為2017年1月—2020年12月，該時間段内SARAL衛星處于漂移軌道，是以資料較為密集。反演使用資料的區域範圍為0°N—20°N，105°E—125°E，輸出重力産品的區域範圍為111.4°N—112.9°N，10°E—12.4°E，在該區域有高精度船載重力、水深測量資料(來自自然資源部)，分辨率達到1′×1′。為了分析SARAL衛星資料的海面分布密度，首先将官方釋出的1 Hz海面高資料在重力場計算的局部海域進行畫圖顯示，分别如圖 2、圖 3所示。
| |
| 圖 2 SARAL衛星1 Hz海面高資料在1′×1′格網内的分布(藍色圓點代表測量值)Fig. 2 Distribution in 1′×1′ grid of 1 Hz sea surface height data from the SARAL satellite (blue dots represent measurements) |
| 圖選項 |
| |
| 圖 3 SARAL衛星1 Hz海面高資料在3′×3′格網内的分布(藍色圓點代表測量值)Fig. 3 Distribution in 3′×3′ grid of 1 Hz sea surface height data from the SARAL satellite (blue dots represent measurements) |
| 圖選項 |
由圖 2、圖 3可知，SARAL衛星1 Hz海面高資料在1′×1′格網内的分布很不均勻，不能完全達到1′×1′分辨率要求，總體上測量點格網分辨率優于3′×3′。實際重力場反演分别按照1′×1′、3′×3′格網大小進行計算，這樣做的原因主要考慮以下3個原因：一是為了與已有1分格網船載重力資料進行比較，友善給出SARAL衛星反演重力場的精度；二是同時為了友善利用國際已經釋出的1分格網的ETOPO-1、DTU18模型；三是考慮到本文方法需要在小區域利用已知格網觀測量進行求解，而理論上格網劃分越細，效果越好，同時考慮到目前試驗區已知的船載水深格網模型為1′×1′格網，是以，采用1′×1′格網進行計算有利于對方法進行驗證。計算中為避免遠區項影響，采用移去恢複方法，先驗模型采用2160階EGM2008模型，按照式(1)—式(3)計算得到南海區域擾動重力如圖 4所示。
| |
| 圖 4 SARAL衛星測高資料反演獲得的1′格網擾動重力資料Fig. 4 The 1′ grid disturbing gravity data obtained from SARAL satellite altimetry data |
| 圖選項 |
經與該區域船載重力資料比較，利用SARAL衛星測高資料及EGM2008模型得到的1′×1′格網擾動重力資料精度達到5.5 mGal，3′×3′格網擾動重力資料精度達到4.8 mGal(表 1)，符合現有衛星測高反演重力場的正常水準。
表 1 SARAL衛星測高資料反演擾動重力精度Tab. 1 Estimation accuracy of disturbing gravity inversion from SARAL satellite
3.2 區域海底地形反演
試驗區域和重力資料反演範圍一緻，具體為10°N—12.4°N，111.4°E—112.9°E，試驗區平均海深3828 m，最大深度4316 m，最小深度1549 m，變化幅度562 m。采用的重力資料為上一節得到的1′×1′格網擾動重力資料。先驗海底地形模型分别采用美國國家環境資訊中心釋出的1′×1′格網ETOPO-1水深模型(與船載水深比較結果見表 2)和丹麥科技大學釋出的1′×1′格網DTU18海深模型(與船載水深比較結果見表 2)，船載水深檢核格網資料在自然資源部2018—2019年的南海實測資料(測線間距約2 km)基礎上格網化得到，分辨率達到1′。
表 2 先驗水深模型的精度統計Tab. 2 Estimation accuracy of seabed topography based on prior bathymetric model
由于曲面函數需要估計5個參數，這就使得重力觀測格網資料有最小解算量，過少的觀測量估計的參數準确度偏低，過多的觀測量又會降低參數的适應範圍，進而降低對高頻資訊的恢複。為了考察不同格網大小劃分對計算結果的影響，将區域格網劃分設定為3×3格網(即每個小區域有9個1′×1′格網)、4×4格網、5×5格網、6×6格網、9×9格網、10×10網格和12×12網格，不同格網劃分條件下的計算結果見表 3，其中基于10×10網格劃分反演得到的海底地形三維圖形如圖 5所示。
表 3 基于ETOPO-1先驗水深模型的海底地形估計精度Tab. 3 Estimation accuracy of seabed topography based on ETOPO-1 prior bathymetric model
| |
| 圖 5 基于ETOPO-1先驗水深模型估計得到的海底地形Fig. 5 Seabed topography estimated based on ETOPO-1 prior bathymetric model |
| 圖選項 |
基于上述方法及DTU18先驗水深模型得到的估計海底地形，與船載資料比較得到結果見表 4，其中基于9×9網格劃分反演得到的海底地形三維圖形如圖 6所示。
表 4 基于DTU18先驗水深模型的海底地形估計精度Tab. 4 Estimation accuracy of seabed topography based on DTU18 prior bathymetric model
| |
| 圖 6 基于DTU18先驗水深模型估計得到的海底地形Fig. 6 Seabed topography estimated based on DTU18 prior bathymetric model |
| 圖選項 |
從以上結果可以看出，雖然理論上大于5個觀測量的重力資料都可以進行高斯曲面估計，但實際上，受制于重力資料本身分辨率及高斯曲面函數的非線性特征，其5個參數的求解受到資料量及資料分布影響還是比較大的，在網格資料劃分上，數量過少或過多都不能取得好的效果。此外，如果海底地形出現複雜的變化趨勢如非光滑、非連續界面，則高斯曲面函數受限于重力觀測資料的分辨率難以重構出複雜的海底地形特征。綜合分析看，對于ETOPO-1先驗模型，10×10格網劃分取得的效果最好，估計得到的海底地形精度相對于ETOPO-1提高了約10 m。對于DTU18模型，9×9格網劃分取得的效果最好，估計得到的海底地形精度相對于DTU18提高了約9 m。綜合來看，對于1′×1′格網海底地形反演而言，按照9×9格網劃分進行分區求解，在不依賴外部海深資料支援下可以取得較優結果。
4 結語
本文從重力資料和海深資料的随機相關特征出發，探索了利用衛星測高反演重力資料及曲面特征函數估計海底地形的方法，得到以下有益結論。
(1) SARAL衛星測高1 Hz資料在南海區域的總體格網密度優于3′×3′，但不能滿足1′×1′格網分辨率的要求，利用船載重力資料進行評估，其反演3′×3′格網重力場精度達到4.8 mGal，反演的1′×1′格網重力場精度優于5.5 mGal。
(2) ETOPO-1、DTU18都是國際上釋出的比較權威的海深模型，其模型建構本身就采用了衛星反演水深和船載測量水深等多種來源的資料，以上述模型為先驗模型，利用衛星測高獲得的重力資料在9×9格網劃分形式下通過高斯曲面函數獲得特征參數進而估計海底地形，其精度較先驗海深模型有約10 m量級的提升，這說明利用最新重力觀測資料在不依賴外部實測水深資料條件下仍然可以對ETOPO-1、DTU18這些已有模型進行精化，也說明利用高斯曲面函數特征參數進行估計的思路是可行的。
(3) 總體而言，本文方法更适合在那些船載測量資料無法有效覆寫或覆寫不均勻的海域，在這些區域可利用高品質的重力資料并結合先驗資訊進行優化提升。實際計算中，高斯曲面函數的估計精度其實和觀測資料的密集程度及精度有關，如果觀測資料分辨率很高，則高斯函數的适應性也會很強，曲面函數就越能反映海底地形區域變化的細部特征(高頻資訊)。是以，重力場越精細，理論上所能反映海底地形的細部特征就越明顯。對于未來以海洋測繪為任務目标的衛星測高技術發展而言，在提高精度的同時更應提高海面高的真實實體分辨率，這樣才有可能利用高分辨率重力場去提升未測海域精細海底地形的反演能力。
作者簡介
第一作者簡介：翟振和(1980—), 男, 博士, 副研究員, 研究方向為空間大地測量與垂直基準。E-mail: [email protected]
初審：張 琳複審：宋啟凡
終審：金 君
資訊