馮進凱，王慶賓，趙東明，楊 洋，黃子炎，黃 炎

（1. 信息工程大學，鄭州 450001；2. 鄭州大學，鄭州 450001）

重力數(shù)據是研究重力場理論和應用的基礎。為使用的方便，重力數(shù)據一般需要以等經緯度格網的數(shù)字形式存儲，以格網中點重力值表征網格區(qū)域內重力場特征。但在實際測量過程中，實測點通常按照測線方式推進、閉合。而且測量過程中容易受到地形等外界因素的干擾，實際測量數(shù)據不一定是格網中點的重力值，這就需要對實測重力點進行格網化，使其按照一定的分辨率等經緯度排列，以利于重力數(shù)據的建模和應用[1]。

傳統(tǒng)的重力數(shù)據格網化的本質是數(shù)學擬合，其原理是對已知的離散重力實測數(shù)據進行數(shù)學內插或外推，從而求得未知點（格網點中點）的重力值。常用的方法有曲面擬合法、反距離加權法[2]、Shepard 擬合方法[3]、Kriging 方法[4,5]以及最近興起的BP-神經網絡和深度學習擬合算法[6]，這些方法從二維數(shù)值擬合的角度出發(fā)處理重力實測數(shù)據，達到了一定的精度，但是這類方法忽略了重力場數(shù)據本身具有的物理特性。文獻[7]研究表明，由于地球重力場信號的中高頻部分和地形之間存在較強的函數(shù)關系，在一定范圍內，地形數(shù)據和重力數(shù)據展現(xiàn)出極強的相似性。而且，現(xiàn)階段全球地形數(shù)據已經能達到相當?shù)木群头直媛?，這為利用結合地形數(shù)據進行重力數(shù)據的格網化奠定了基礎?；诖耍疚睦弥亓鰯?shù)據頻譜疊加原理，在顧及重力數(shù)據物理特性的基礎上結合實驗區(qū)地形數(shù)據，提出一種基于“移去-恢復”理論的重力數(shù)據格網化方法，設計出“計算-移去-推估-恢復”的四步重力數(shù)據格網化方法（簡稱“移去-恢復”法），實驗結果表明，相比于傳統(tǒng)的Kriging 格網化方法，本文提出的方法在精度上有一定的提升，算法穩(wěn)定性更高，自洽性更好，而且格網化結果能夠反映目標區(qū)域更多的局部特征。

1 原理和方法

1.1 均衡理論

根據均衡理論[8]，均衡異常應滿足浮平衡理論，即均衡重力異常ΔgI理論上應滿足如下條件：

其中， Δ g空為空間重力異常，δ gB和 δIC分別為布格改正和地形均衡改正。布格改正δ gB可表示為：

其中，h 為高程， - 0 .1116h 為層間改正，δgTC為局部地形改正。在局部平面坐標系下，局部地形改正計算時，為了兼顧計算精度和效率，通常將積分區(qū)域σ 劃分為中央區(qū) σ0和遠區(qū)σ -σ0，中央區(qū)采用棱柱積分法，遠區(qū)采用線性卷積方式逼近真實積分[9]：

其中，( x , y , z )為流動點平面坐標，G 為萬有引力常量，δ 為地球密度常數(shù)，h 和 hp分別為流動點和計算點高程。

均衡改正采用較為符合地球實際均衡狀態(tài)的Airy—Heiskanen 模型[10,11]，該理論認為常密度的地殼漂浮在密度為 ρm的地幔上并保持平衡，超出海平面的部分或低于海平面的部分分別在均衡面以下的部分以“山根”或“反山根”的形式進行補償。在局部平面直角坐標系下，直接給出均衡改正的計算表達式：

其中，z2=- T，z1=-（ T+t） ，T 為均衡面深度，t 為補償深度， Δδ 為地殼和地幔部分的密度差，均衡改正計算方法參照局部地形改正。一般來說，空間重力異常變化較為劇烈，在構造高分辨的格網數(shù)據時對測點的密度要求較高，而均衡重力異常變化較為平緩，對測點的密度要求較低，比較適合重力點的格網化。

1.2 Kriging 法

Kriging 方法是一種常用的擬合方法，廣泛的應用于數(shù)據處理以及工程應用等領域，從本質上講，該方法是利用區(qū)域原始數(shù)據對估計點進行的最優(yōu)無偏估計[4,5]。其滿足以下條件：

其中，E 與C 代表數(shù)學期望與方差，N 與 Ni分別表示計算點與擬合點的大地水準面高。插值模型如下：

其中， λi為每個擬合點的權值。通常利用變異函數(shù)代替方差元素：

其中，h 表示滯后距，Num(h) 表示滯后距個數(shù)，詳細推導見文獻[5]。

1.3 “移去-恢復”方法

“移去-恢復”理論在局部重力場建模中有著廣泛應用[12,13]，基本思想為：根據重力場頻譜疊加性原理將重力場信號分為低頻、中頻、高頻三部分，重力場元看作不同頻段重力場信號的組合。隨著重力衛(wèi)星技術的發(fā)展，低頻信號已能達到相當精度，高頻和超高頻部分主要來自地形起伏影響，可以通過高分辨率的地形數(shù)據計算得到。在計算時，移去輸入數(shù)據的低頻、高頻信號，對殘余值進行建模計算，最后在計算結果中恢復相應頻段的重力場信號，即可得到計算點完整頻段的重力場元數(shù)據[4,14]，其計算步驟可簡單概括如下：

（1）首先在重力場元中移除低頻和高頻頻段信息，獲得殘差重力場元；

其中，L表示擾動位T 的泛函，TRes為殘余擾動位，TM和TT分別為低頻信號和高頻信號。

（2）利用步驟（1）中產生的殘差重力場信息結合適當?shù)慕７椒ㄟM行模型構建；

（3）最后，在計算點中恢復步驟（1）中移去的低頻和超高頻信息。

需要注意的是，“移去-恢復”過程中，移去和恢復的部分必須是可以在一定的精度條件下精確計算得到的，否則“移去-恢復”過程將無法完成。

同理，公式(1)可改寫為如下形式：

空間重力異?？梢钥醋魇遣几窀恼?、均衡改正和均衡異常三種信號的疊加，從1.1 中可以看出，布格改正和均衡改正均與地形數(shù)據相關，可通過地形數(shù)據和地殼、地幔密度信息計算得到，符合“移去-恢復”理論的基本前提?，F(xiàn)將“移去-恢復”理論應用到重力數(shù)據格網化中，基本思路可概括為“計算-移去-推估-恢復”四步，具體計算步驟如下：

（1）計算：根據Airy 均衡理論，利用式(2)-(4)并結合高精度、高分辨率的地形數(shù)據和密度信息計算目標區(qū)域內布格改正和均衡改正的改正項；

（2）移去：利用式(1)在建模點的空間重力異常中移去地形數(shù)據產生的影響（布格改正、均衡改正），得到在建模點較為平滑的均衡重力異常；

（3）推估：利用1.2 中介紹的Kriging 格網化方法對均衡重力異常進行格網化，獲取目標點的均衡重力異常數(shù)據；

（4）恢復：利用步驟（1）計算的數(shù)據，在目標點均衡重力異常的基礎上恢復地形數(shù)據的影響，即可得到目標點的空間重力異常，公式如下：

2 數(shù)值實驗

實驗區(qū)位于南非共和國境內的馬塔貝萊高原，范圍如圖1 所示，實驗區(qū)內地形變化劇烈，重力場信號復雜，地形數(shù)據采用美國國家地球物理數(shù)據中心(NGDC)開發(fā)的Etopo1 的1′分辨率格網化地形數(shù)據[15]，實驗區(qū)內最大高程落差 2334 m，平均高程1052.2284 m，選用離散點所在地形格網為離散點高程，顧及地形改正和均衡改正邊界效應，選取圖1 紅框部分為中心計算區(qū)域。從“GRAVCD-africa”數(shù)據集中搜集到計算區(qū)1789 個實測空間重力異常數(shù)據（圖1 中三角點和圓點），地形數(shù)據和實測點空間重力異常數(shù)據統(tǒng)計信息見表1。

圖1 實驗區(qū)域高程及點位分布示意圖Fig.1 Elevation and points distribution of the experimental area

表1 高程和實測點空間重力異常統(tǒng)計表Tab.1 Statistics table of free-air gravity anomaly and elevation

為驗證本文所提方法的正確性，隨機選擇計算區(qū)內894 點為建模點（圓點），其余895 個點為檢核點（三角點），設計如下對比實驗：

方案一：直接格網化方法，利用傳統(tǒng)的格網化方法直接對實驗區(qū)內的建模點進行數(shù)據格網化，格網化時采用物理大地測量學中常用的 Kriging 插值方法[17,18]，其為區(qū)域原始數(shù)據對估計點進行的最優(yōu)無偏估計。利用離散建模點位置信息和空間重力異常得到檢核點空間重力異常數(shù)據，最后和檢核點實測數(shù)據進行對比，從而得到本方案的格網化精度；

方案二：利用本文提出的“計算-移去-擬合-恢復”格網化方法，并結合實驗區(qū)高精度、高分辨率的地形數(shù)據和密度信息構造計算區(qū)域的均衡改正、地形改正和層間改正項，在建模點移去各項改正，得到建模點均衡重力異常，利用Kriging 格網化方法將實驗區(qū)格網化均衡異常格網化；最后根據點所在網格恢復移去項（均衡改正和布格改正），得到檢核點空間重力異常，并和檢核點實測數(shù)據進行對比，分析格網化精度。

2.1 方案二各類改正項分析

方案二“計算”步驟中，均衡面深度選為40 km，地殼密度為2670 kg/m3，巖漿密度為3270 kg/m3，地形改正積分半徑選取30′，利用以上數(shù)據分別計算出中心區(qū)域地形改正、均衡改正和層間改正，計算結果見圖2、表2。

圖2 計算區(qū)各項改正示意圖Fig.2 Diagram of various corrections

表2 各項改正統(tǒng)計表/mGalTab.2 Statistical table of all types of corrections/mGal

分析圖2 可以發(fā)現(xiàn)，在計算區(qū)內地形改正、均衡改正和層間改正的分布和地形趨勢基本一致，相似度極高，而且計算區(qū)內任何位置的改正項均可以通過嚴密公式和高分辨率的地形數(shù)據以及地殼密度項計算出來，符合“移去-恢復”理論的基本條件。

方案二“移去”過程中，依次移除建模點空間重力異常中的布格改正、均衡改正，分別得到布格重力異常和均衡重力異常，統(tǒng)計信息見表3，繪制如下，見趨勢圖3（為方便對比，圖中數(shù)據消除均值）。

圖3 “移去”過程趨勢圖Fig.3 Diagram of “Remove” process

表3 建模點各項重力異常統(tǒng)計表/mGalTab.3 Statistics table of all types of gravity anomaly/mGal

從圖3 中可以看出，在空間重力異常中依次加入布格改正和均衡改正后，建模點的重力異常信號變化由原來的23.93 mGal 衰減為12.8662 mGal，差值閾由原來的140.3415 mGal 減小為78.8301 mGal，毛刺點明顯減少，變化趨勢更加平滑，所以空間重力異常在經過布格改正和均衡改正后更適合目標點的格網化，數(shù)據處理過程中的誤差將會減小。

2.2 不同格網化方法精度分析

按照方案一和方案二組織對比實驗，統(tǒng)計兩個實驗方案的結果于表4，差值見圖4。

圖4 兩種方法各點差值示意圖Fig.4 Diagram of the differences between the two methods

表4 兩種方法差值統(tǒng)計表/mGalTab.4 Statistical table of differences between two methods/mGal

實驗結果證明，無論哪種方案均能以一定的精度推估未知點的空間重力異常，其中方案一從數(shù)學角度建立已知點和未知點的聯(lián)系，進而根據建模點的空間重力異常推估出檢核點空間重力異常，檢核精度為5.5168 mGal；方案二在數(shù)學格網化的基礎上顧及了地球自身的物理信息，擬合精度為3.9996 mGal，相比于方案一，精度提升了27.5%，從精度上來看，本文提出的方法更優(yōu)。“移去-恢復”法計算結果差值最大值和最小值相比傳統(tǒng)方法均有了一定的收斂，閾值縮小近一倍。從圖4 中看，方案二各檢核點差值分布更加集中，誤差分布更加均勻，粗差點更少。這是由于“移去-恢復”中在純數(shù)學擬合中加入了地球的物理信息，整體擬合趨勢和重力數(shù)據分布趨勢更加符合，擬合精度更好，差值分布更收斂，殘差值波動也更小。

根據以上兩種實驗方案，利用計算區(qū)內全部的1789 個重力實測數(shù)據作為建模數(shù)據，構造出計算區(qū)內1′分辨率的空間重力異常，得到實驗結果如圖5 所示，再將1789 個建模數(shù)據作為檢核數(shù)據驗證算法的自洽性，統(tǒng)計結果見表5 和圖6。

圖5 兩種方法構造的區(qū)域重力異常圖Fig.5 Regional Free-air gravity anomaly constructed by two methods

圖6 和表5 表明，本文提出的“計算-移去-推估-恢復”方法自洽性檢核的各項數(shù)學統(tǒng)計指標均優(yōu)于直接推估方法，精度更高，系統(tǒng)差更小。而且在圖5 中可以看出，“移去-恢復”法所構造的局部空間重力異常包含更多的細節(jié)特征，細部信息刻畫更為詳細。這是 因為傳統(tǒng)的數(shù)學推估方法在計算時，會損失數(shù)據中原有的高頻信號，推估結果較為平滑；本文所提方法在“移去-恢復”的過程中，充分考慮了地形數(shù)據的影響，補齊了實驗區(qū)內空間重力異常的高頻項，使得重力場元的頻段更加完整。

圖6 兩種方法自洽性差值示意圖Fig.6 Diagram of self-consistent error by two methods

表5 兩種方法自洽性誤差值統(tǒng)計表/mGalTab.5 Statistics table of self-consistent error by two methods/mGal

3 結 論

本文將“移去-恢復”理論應用到局部重力數(shù)據格網化中，在格網化過程中加入地形數(shù)據和密度信息，將傳統(tǒng)的純數(shù)學擬合轉化為顧及地球物理場信息的格網化，通過實驗區(qū)實測數(shù)據對比實驗，得到如下結論：

（1）實測空間重力異常較為“毛糙”，直接對其進行數(shù)學擬合時，過程性誤差較大；在消除掉布格改正和均衡改正后數(shù)據變得平滑，有利于重力數(shù)據的格網化；

（2）實驗區(qū)內，“移去-恢復”方法構造的目標區(qū)域空間重力異常精度相比于Kriging 格網化的精度更高，差值波動范圍更小，而且能夠描述實驗區(qū)內重力場的細節(jié)特征，重力場元的頻段更加完整；

（3）“移去-恢復”方法利用地形數(shù)據將變化劇烈的空間重力異常轉化為較為平滑的均衡重力異常，然后再進行數(shù)據處理?？梢灶A測，當目標區(qū)域內實測重力數(shù)據較少時，“移去-恢復”結合高分辨率地形數(shù)據方法的優(yōu)勢將更加顯著。