一種基于“移去-恢復(fù)”理論的重力數(shù)據(jù)格網(wǎng)化方法

馮進(jìn)凱，王慶賓，趙東明，楊 洋，黃子炎，黃 炎

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

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

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

1 原理和方法

1.1 均衡理論

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

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

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

其中，( x , y , z )為流動(dòng)點(diǎn)平面坐標(biāo)，G 為萬有引力常量，δ 為地球密度常數(shù)，h 和 hp分別為流動(dòng)點(diǎn)和計(jì)算點(diǎn)高程。

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

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

1.2 Kriging 法

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

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

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

其中，h 表示滯后距，Num(h) 表示滯后距個(gè)數(shù)，詳細(xì)推導(dǎo)見文獻(xiàn)[5]。

1.3 “移去-恢復(fù)”方法

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

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

其中，L表示擾動(dòng)位T 的泛函，TRes為殘余擾動(dòng)位，TM和TT分別為低頻信號(hào)和高頻信號(hào)。

（2）利用步驟（1）中產(chǎn)生的殘差重力場(chǎng)信息結(jié)合適當(dāng)?shù)慕７椒ㄟM(jìn)行模型構(gòu)建；

（3）最后，在計(jì)算點(diǎn)中恢復(fù)步驟（1）中移去的低頻和超高頻信息。

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

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

空間重力異?？梢钥醋魇遣几窀恼⒕飧恼途猱惓ＨN信號(hào)的疊加，從1.1 中可以看出，布格改正和均衡改正均與地形數(shù)據(jù)相關(guān)，可通過地形數(shù)據(jù)和地殼、地幔密度信息計(jì)算得到，符合“移去-恢復(fù)”理論的基本前提。現(xiàn)將“移去-恢復(fù)”理論應(yīng)用到重力數(shù)據(jù)格網(wǎng)化中，基本思路可概括為“計(jì)算-移去-推估-恢復(fù)”四步，具體計(jì)算步驟如下：

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

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

（3）推估：利用1.2 中介紹的Kriging 格網(wǎng)化方法對(duì)均衡重力異常進(jìn)行格網(wǎng)化，獲取目標(biāo)點(diǎn)的均衡重力異常數(shù)據(jù)；

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

2 數(shù)值實(shí)驗(yàn)

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

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

表1 高程和實(shí)測(cè)點(diǎn)空間重力異常統(tǒng)計(jì)表Tab.1 Statistics table of free-air gravity anomaly and elevation

為驗(yàn)證本文所提方法的正確性，隨機(jī)選擇計(jì)算區(qū)內(nèi)894 點(diǎn)為建模點(diǎn)（圓點(diǎn)），其余895 個(gè)點(diǎn)為檢核點(diǎn)（三角點(diǎn)），設(shè)計(jì)如下對(duì)比實(shí)驗(yàn)：

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

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

2.1 方案二各類改正項(xiàng)分析

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

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

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

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

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

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

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

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

2.2 不同格網(wǎng)化方法精度分析

按照方案一和方案二組織對(duì)比實(shí)驗(yàn)，統(tǒng)計(jì)兩個(gè)實(shí)驗(yàn)方案的結(jié)果于表4，差值見圖4。

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

（2）實(shí)驗(yàn)區(qū)內(nèi)，“移去-恢復(fù)”方法構(gòu)造的目標(biāo)區(qū)域空間重力異常精度相比于Kriging 格網(wǎng)化的精度更高，差值波動(dòng)范圍更小，而且能夠描述實(shí)驗(yàn)區(qū)內(nèi)重力場(chǎng)的細(xì)節(jié)特征，重力場(chǎng)元的頻段更加完整；

（3）“移去-恢復(fù)”方法利用地形數(shù)據(jù)將變化劇烈的空間重力異常轉(zhuǎn)化為較為平滑的均衡重力異常，然后再進(jìn)行數(shù)據(jù)處理。可以預(yù)測(cè)，當(dāng)目標(biāo)區(qū)域內(nèi)實(shí)測(cè)重力數(shù)據(jù)較少時(shí)，“移去-恢復(fù)”結(jié)合高分辨率地形數(shù)據(jù)方法的優(yōu)勢(shì)將更加顯著。