海域航空重力快速構建區域大地水準面

王晨陽， 駱遙， 熊盛青， 劉詩華， 劉祖鑒

1 中國地質大學(北京) 地球物理與信息技術學院， 北京 100083 2 中國自然資源航空物探遙感中心， 北京 100083 3 國防科技大學智能科學學院， 長沙 410073

0 引言

航空重力以航空器為測量平臺，綜合利用重力傳感器、慣性導航系統、全球衛星定位系統等測定空中重力加速度(孫中苗等，2004；熊盛青等，2010).航空重力測量始于20世紀 80年代(張永明等，2006)，隨著導航技術以及高精度重力傳感器技術的應用，近年來航空重力測量技術獲得了突破性進展，測量精度達1×10-5～2×10-5m·s-2(熊盛青，2009)，空間分辨率3～5 km.航空重力測量不受作業區域限制，能夠快速獲得均勻分布的高精度重力場觀測資料，同船測、地面重力測量相比作業效率更高.海域航空重力測量可以彌補衛星測高數據在近海及海陸交互帶、浮冰區等目標區域觀測精度較低的問題(Brozena et al.，1990；Hwang et al.，2006)，是海洋地球物理探測的重要手段之一.

大地水準面是與全球無潮平均海平面重合并延伸至大陸內部的重力等位面，是空間高程的基準面(Grafarend，1994；吳曉平，2006；魏子卿，2009；林淼等，2012)，其形狀反映了地球內部物質密度結構、應力場、地球旋轉和洋流等信息(Vanichek and Christou，1994).大地水準面與地球重力場模型關系密切，現代大地水準面的確定以地球重力場的位理論為基礎(張利明和李斐，2005)，近年來通過航空重力測量精化大地水準面成為研究熱點(蔣濤，2013；孫中苗等，2014；劉金釗等，2020).目前，我國已在管轄海域開展了大規模的航空重力勘探，為海域基礎地質調查和油氣資源調查提供了重要基礎性資料(李文勇等，2010；張玄杰等，2016).利用已有航空重力數據構建管轄海域大地水準面，對拓展航空重力勘探應用方向具有現實意義.本文針對海域航空重力勘探實際，利用高精度向下延拓等關鍵技術，實現了一種快速高效構建海域大地水準面的方法，可為確定海域大地水準面提供技術支撐.

1 大地水準面建模

大地水準面與重力場密切相關，Brun公式(Heiskanen and Moritz，1967)給出了二者的直接聯系，即大地水準面起伏可以由重力擾動位除以正常重力值確定.大地水準面建模的關鍵在于將實測航空重力異常轉換為重力擾動位.

1.1 位場轉換近似

大地測量中通常通過斯托克斯積分計算水準面(許厚澤，2006)，即

(1)

其中：N是地心緯度為φ、地心經度為λ處的水準面起伏，R是地球半徑，S(ψ)是斯托克斯積分核，Δg是球面積分區域σ對應的重力異常.如果該積分區域有限，可將球面近似為平面，即斯托克斯積分將近似為

(2)

其中：N是(x，y)處的水準面起伏，γ是正常重力值，Δg是平面積分區域s上相應的重力異常.該積分是卷積形式，利用傅里葉變換可將其寫成頻率域形式，即

(3)

其中，F[·]表示傅里葉變換，u、v是Δg傅里葉變換后x、y方向對應的波數，F[·]-1表示傅里葉變換的逆變換.因此，大地水準面起伏N可以通過公式(3)表達的頻率域位場轉換確定.

對于航空重力勘探，測量高度通常高于積分面，需要將測量的重力異常向下延拓至積分面.陸域航空重力測量觀測面多是隨地形起伏的曲面，難以嚴格應用位場轉換條件進行向下延拓.海域航空重力通常測量的是距海平面或橢球面等高的重力異常，具備平面位場轉換的適應條件，可直接向下延拓至積分面，即

(4)

其中：Δgh是實測航空重力異常(空間重力異常)，h>0是測量高度，u、v是Δgh傅里葉變換后x、y方向對應的波數.因此，海域水準面起伏可利用航空重力資料通過重力勘探中位場轉換方法近似計算，即對航空重力進行向下延拓，再對其進行垂向積分并除以正常重力值獲得.由于位場轉換要求待轉換的位場為平面網格數據，利用公式(3)和公式(4)確定大地水準面在一定范圍內忽略了地球曲率影響，但頻率域位場轉換算法處理效率非常高，可用于航空重力測量成果快速構建局部海域大地水準面.

1.2 高精度向下延拓技術

海域航空重力采用固定飛行高度測量，可直接將航空重力異常從平均飛行高度(橢球高)向下延拓至橢球面.向下延拓可以從數學上拉近觀測平面與場源間的距離，改善信噪比，增強地質信息，但重力異常向下延拓卻是經典的不適定問題(Blakely et al.，1996；曾華霖，2005；黃謨濤等，2018)，采用公式(4)進行向下延拓將強烈的表現出對數據高頻部分的振蕩，向下延拓后的數據通常無法直接使用.由于還需對向下延拓的重力異常進行垂向積分，也可將頻率域向下延拓算子同垂向積分算子連乘進行位場轉換，該算子為

(5)

由于公式(5)的位場轉換算子隱含垂向積分操作，相比直接向下延拓，一定程度上會減少對異常高頻部分放大的程度，略顯穩定.Cooper(2004)曾應用類似(5)式的位場轉換算子進行處理，再進行空間域求導得到相對穩定的向下延拓結果.盡管公式(5)表示的位場轉換較直接向下延拓更穩定，但直接對實測航空重力數據進行處理，仍將出現高頻振蕩，特別是航測高度較高的情況下.目前，國內海域航空重力勘探測量高度為800 m或1000 m，通常的向下延拓處理仍難以獲得橢球面上穩定的異常數據.

位場轉換還存在假值填補與數據擴邊問題.快速傅里葉變換(FFT)對網格數據的行、列數有一定要求，通常要求為2的自然數次冪.網格化后的航空重力異常數據一般不滿足該條件，需要對網格數據進行擴邊，擴邊效果直接影響后續處理的精度，Nagy和Fury(1990)曾就數據擴邊對重力異常轉換大地水準面精度的影響進行過討論.實際航空重力異常資料還因測網設計、偏離航線、測線掐線等限制，網格化后某些區域無法形成相應的內插數據進而存在網格格點空白，處理時要求對空白區進行填補，這也將影響處理的精度.為了解決上述問題并提高處理精度，本文采用位場向下延拓的最小曲率方法對航空重力異常進行向下延拓(駱遙和吳美平，2016).具體采用如下迭代形式進行向下延拓，即

(6)

其中：Δgk表示第k次迭代后的航空重力向下延拓結果；η是權重系數；downcα(h)表示正則化參數為α的最小曲率向下延拓算子，h>0是向下延拓距離；upc(h)是延拓距離為h的向上延拓算子.上述處理解決了實際非規則航空重力測量中空白區及網格擴邊問題，提升了航空重力勘探向下延拓處理的實用性，可高精度的獲取下延平面上觀測異常對應的位函數.

1.3 移去與恢復

上述向下延拓獲得的異常即為擴邊后的下延重力異常，直接對其積分即可得到相應的位函數.對異常進行積分時，零頻u=v=0處算子存在奇異，通常將奇異時算子設為0(Lourenco et al., 1973)，加之存在積分誤差，位函數的長波長部分無法準確確定.因此，位場轉換所確定的位函數是相對的，而擾動位的基準直接關系到水準面，必須確定其起算基準.為解決該問題，這里采用移去-恢復法的思想(黃志洲等，2004)，用EGM2008等重力場模型確定低頻重力場，其最終確定的大地水準面包括局部航空重力異常計算部分和重力場模型計算部分，即

N=NΔg+NGM，

(7)

其中：N是大地水準面起伏，NΔg是局部航空重力異常確定的水準面起伏，NGM是重力場模型確定的大地水準面.

2 模擬計算

為了檢驗構建大地水準面的位場轉換方法，我們采用EMG2008模型(階次2159)模擬航空重力數據和大地水準面，檢驗方法的精度.按1∶1000000標準圖幅范圍，計算我國香港幅(114°E—120°E，20°N—24°N)海域高1 km(大地高)航空重力測量數據，圖1給出了模擬的海域航空重力測量數據，其中網格距為1 km，對應陸域范圍內的數據均設為空白.航空重力測量通過引入重力基點可以獲取絕對重力場，圖1a的海域航空重力場數據為EGM2008計算的絕對重力場.我們將EGM96模型(階次360)對應的絕對重力場移去，圖1b給出了移去后的航空重力異常.采用1.2節向下延拓方法將圖1b重力異常向下延拓1 km歸算至橢球面，并對延拓后的重力異常進行積分換算成位函數，計算相應的水準面起伏(圖2a).圖2b給出了通過(7)式最終確定的大地水準面，其中模型水準面部分采用EGM96.我們將計算的大地水準面同EGM2008模型進行比較，二者最大相差37 cm，標準差10.27 cm， EGM2008模型全球平均精度為13 cm(Pavlis et al.，2012)，構建的水準面精度與此相當，說明該方法具有較高的計算精度.

圖1 模擬的海域航空重力場(a)及移去后航空重力異常(b)Fig.1 The simulated airborne gravity field of sea area (a) and airborne gravity anomaly after removal (b)

3 實際航空重力資料計算

為檢驗位場轉換方法實際構建海域大地水準面效果，我們采用航空重力實測資料進行處理.使用的數據為美國國家海洋和大氣管理局(National Oceanic and Atmospheric Administration，NOAA)所屬的國家大地測量局(National Geodetic Survey，NGS)在波多黎各附近加勒比海海域實測航空重力數據.該航空重力測量是美國構建下一代垂直基準GRAV-D計劃的一部分，并利用地面絕對重力測量點將空中重力異常數據進行了歸算，最終航空重力測量成果以測量高度上的絕對重力場表達.圖3給出了實測的航空重力測線及航空重力數據.采用塞斯納獎狀噴氣飛機(Cessna Citation II )搭載美國Micro-g LaCoste公司的TAGS航空重力系統進行測量，設計的飛行高度為10.668 km，實際平均飛行高度為11226.11 m，測線間距為10 km，切割線間距約40 km.圖3b航空重力數據的網格距為2.0 km(航空地球物理測量中網格距通常為測線間距的1/5～1/4).由于測線邊緣不整齊，網格化后的格網數據南北部邊緣存在數據空白.此外，區內兩條測線上部分存在問題的數據已被掐去，造成網格化后的格網數據開“天窗”.圖3反映了實際航空重力測量的典型特征，目前一些通過模擬航空重力數據精化水準面的研究中通?；乇芰松鲜鰯祿卣?，直接用EGM2008等重力場模型來模擬不含數據空白的規則矩形格網數據.根據1.3節的技術路線，利用EGM2008重力場模型(階次2159)計算重力場的模型部分，并將其從實測數據中移去，圖4a給出了移去后的重力異常，異常幅度為-8.3×10-5～7.87×10-5m·s-2.利用1.2節的方法將圖4a的重力異常向下延拓，圖4b給出了下延后的重力異常.向下延拓的距離為11226.11 m(實測平均橢球高)，圖4b可視為圖4a歸算至WGS-84橢球面上的重力異常.圖4b經向上延拓至原測量面并同圖4a比較，誤差的標準差為0.26×10-5m·s-2，99%的網格點誤差幅度不高于0.91×10-5m·s-2，說明向下延拓具有很高的數值精度.

圖3 航空重力測線(a)及實測航空重力數據(b)Fig.3 Airborne gravity survey lines (a) and airborne gravity field data (b)

圖4 航空重力異常(a)及下延至WGS-84橢球面上的異常(b)Fig.4 Airborne gravity anomaly (a) and its downward continuation result on the WGS-84 ellipsoid (b)

將圖4b歸算至WGS-84橢球面上的重力異常轉化為相應的擾動重力位，圖5a給出了相應的計算結果.擾動位精度直接影響大地水準面精度，為了消除積分的水平誤差，將圖5a的位函數反算為航空重力異常，并將其同圖4a進行比較，圖5b給出了位函數轉換重力異常的誤差.圖5b給出的誤差分布表明計算的位具有很高的精度，誤差呈正態分布，誤差幅度約2×10-5m·s-2量級，標準差為0.27×10-5m·s-2，基本同向下延拓誤差的標準差0.26×10-5m·s-2相當.盡管誤差的標準差很小，但誤差的平均值為2.11×10-5m·s-2，表明零頻處理中存在相應系統差，需進一步消除.我們將圖5a的位函數求導，轉換為重力異常同圖4b進行比較，二者誤差的標準差為5×10-12m·s-2，幾乎為0，誤差的平均值為2.11×10-5m·s-2，這表明零頻置0的處理中忽略了2.11×10-5m·s-2的重力場直流分量.因此，零頻處理中引起的系統誤差造成圖5a擾動位比實際要大，其忽略的重力場常量在垂向的積分近似為0.2369 m2·s-2(即2.11×10-5m·s-2×11226.11 m).

補償圖5a資料中0.2369 m2·s-2的系統差后，計算剩余異常所引起水準面起伏，圖6a給出了圖4a對應的水準面起伏，其幅度最大約30 cm.最終，利用公式(7)將圖6a同EGM2008水準面融合，圖6b給出了該測區最終確定的大地水準面.

圖5 由剩余異常計算的位函數(a)及誤差(b)Fig.5 Gravity potential field calculated by residual anomaly (a) and error (b)

圖6表明該區的大地水準面起伏在-70.84 m至-40.61 m之間，精化大地水準面起伏不超過38 cm，從位場轉化誤差分析其精度約為2×10-5m·s-2，與測區航空重力實測精度1.85×10-5m·s-2基本相當，這是從位場轉換方法本身的精度去評價構建大地水準面的精度.為了進一步檢驗方法的精度，這里采用測區內美屬維爾京群島上的21個GPS/水準點進行檢驗(由于向下延拓存在過場源問題，該檢驗方法僅為近似評估)，圖7給出了GPS/水準點分布.將GPS/水準點的水準面同EGM2008模型大地水準面及通過航空重力測量精化的大地水準面進行比較，圖8給出各點位的水準面及大地水準面對比情況.GPS/水準點給出的水準面同EGM2008大地水準面模型及本文給出的大地水準面趨勢均較為一致，但GPS/水準測量采用了VIVD09垂直基準，同大地水準面間存在水平差.圖8表明VIVD09的參考水準要高出EGM2008大地水準面2.5 m的量級，如果忽略二者間的水平差，以GPS/水準測量為標準，EGM2008模型的大地水準面精度可用二者差值的標準差代表，其精度為10.85 cm.本文給出的大地水準面同GPS/水準測量也存在類似水平差，以GPS/水準測量為標準，二者差值的標準差為6.36 cm，較EGM2008大地水準面的精度提高了約4 cm.上述檢驗使用的GPS/水準點均位于島嶼處，檢驗中假設了向下延拓過場源時將參考橢球外部的質量壓縮至橢球面上，存在一定的近似.為此，我們采用美國xGEOID18B水準面模型進行了檢驗.圖6b給出的水準面同xGEOID18B相比，最大相差18.9 cm，二者差值的標準差為5.08 cm，可見本文的海域航空重力測量快速精化大地水準面方法具有較高的精度.

圖6 大地水準面修正量(a)及修正后大地水準面(b)Fig.6 Correction of geoid (a) and final geoid determined (b)

圖7 GPS/水準點分布Fig.7 GPS/level point distribution

圖8 水準面與EGM2008大地水準面及航空重力確定的大地水準面對比Fig.8 Geoid of GPS/level compared with the EGM2008 geoid and the geoid determined by airborne gravity data

4 結論

構建高精度大地水準面將是我國航空重力勘探下一步應用的重要方向，基于移去-恢復法思想，本文采用重力勘探中位場轉換技術，特別是引入向下延拓的最小曲率方法，給出了一種快速構建海域大地水準面方法.方法通過移去EGM2008地球重力場，采用最小曲率法對移去后剩余的航空重力異常進行向下延拓并在頻率域對其積分獲取擾動位，最終通過結合EGM2008水準面模型實現對海域大地水準面精化.理論模擬和實際航空重力測量數據處理表明，盡管該方法在構建海域大地水準面存在位場轉換近似，但具備較高的計算效率和精度，可為利用海域航空重力勘探成果快速構建大地水準面提供技術支撐.航空重力勘探屬于相對重力測量，本方法需要結合地球重力場模型，使用的EGM2008等全球重力場模型在一定程度上限制了構建的海域大地水準面精度，我們希望同測繪學界進一步合作，共同構建我國自主的高精度全球重力場模型，進一步提升我國大地水準面的精度.

致謝本文圖3使用了美國NOAA的數據，在此表示感謝.感謝兩位匿名審稿專家對本文提出的批評和建議.