基于拉格朗日中值定理的航空重力異常向下延拓方法研究

胡 起 徐新禹 趙永奇 丁子凡

1 武漢大學(xué)測(cè)繪學(xué)院，武漢市珞喻路129號(hào)，430079

重力異常向下延拓是物理大地測(cè)量的經(jīng)典問題。長期以來，國內(nèi)外開展了很多研究來解決向下延拓中的不適定問題，其中以求解球外Dirichlet問題的逆Poisson積分法應(yīng)用最為廣泛。迭代法是解算逆泊松積分的常用方法。Martine與Huang[1-2]利用迭代方式將加拿大落基山脈的航空重力值進(jìn)行向下延拓；Kingdon、Ma與Zeng等[3-5]對(duì)迭代法進(jìn)行了更深入的研究。正則化方法也是解決病態(tài)問題的常用方法。王興濤等[6]采用Tikhonov正則化方法對(duì)5′×5′的航空重力異常進(jìn)行向下延拓。蔣濤等[7]則將嶺估計(jì)與廣義嶺估計(jì)應(yīng)用于航空重力向下延拓。趙啟龍等[8]使用半?yún)?shù)估計(jì)與最小二乘配置方法對(duì)臺(tái)灣地區(qū)的航空重力進(jìn)行向下延拓。在頻率域進(jìn)行重力延拓也是常用的方法。將重力異常從空間域利用FFT算法轉(zhuǎn)換到頻率域，在頻率域進(jìn)行向下延拓，最后將延拓后的重力異常從頻率域轉(zhuǎn)換回空間域[9]。寧津生等[10]使用小波分析研究了多尺度邊緣約束的重力場向下延拓。黃謨濤等[11]提出基于向上延拓的向下解析延拓方法，利用向上延拓與垂直梯度進(jìn)行航空重力向下延拓。Zhang也提出利用積分中值定理進(jìn)行航空重力向下延拓的數(shù)值方法[12]。不同于逆Poisson積分法，還有一些方法是利用物理信息進(jìn)行重力場延拓的研究，如直接代表法、虛擬點(diǎn)質(zhì)量法、聯(lián)合重力場模型與地形進(jìn)行重力場向下延拓等。

本文結(jié)合拉格朗日中值定理以及重力異常垂直梯度在垂向上的分布特點(diǎn)，改進(jìn)使用重力異常垂直梯度進(jìn)行重力異常延拓的方法，回避了向下延拓的不適定性。針對(duì)經(jīng)典積分求解重力異常垂直梯度的奇異性，本文使用差分的方式計(jì)算航空重力高度面上的垂直梯度，并在臺(tái)灣地區(qū)進(jìn)行模擬與實(shí)測(cè)航空重力異常向下延拓實(shí)驗(yàn)。

1 基本原理與方法

1.1 利用重力異常垂直梯度進(jìn)行向下延拓

重力異常向下延拓中，最直接的方法是利用垂直梯度進(jìn)行向下延拓。根據(jù)拉格朗日中值定理，有：

(1)

將式(1)變形得到：

(2)

重力異常垂直梯度(Δgr)在垂向上的變化較為緩慢，且呈現(xiàn)出與高度較強(qiáng)的線性關(guān)系。為了說明重力異常垂直梯度在垂向的分布特點(diǎn)，本文在臺(tái)灣地區(qū)隨機(jī)選取100個(gè)點(diǎn)，采用2 160階次的EGM2008重力場模型模擬每個(gè)點(diǎn)從海平面到海平面以上10 000 m每隔100 m的重力異常，點(diǎn)位分布見圖1。圖2為各點(diǎn)位上每隔100 m模擬重力異常垂直梯度與高度的散點(diǎn)圖。

圖1 Δgr垂向分布實(shí)驗(yàn)的點(diǎn)位分布Fig.1 Distribution of points validating the radial pattern of Δgr

圖2 各點(diǎn)位上Δgr與高度的散點(diǎn)圖Fig.2 Scatter point plots of vertical gradients and height at each test point

由圖2可以看出，多數(shù)點(diǎn)位上的模擬重力異常垂直梯度與高度呈現(xiàn)出強(qiáng)線性關(guān)系。為進(jìn)一步說明這種線性關(guān)系，求得各點(diǎn)位上每隔100 m模擬重力異常垂直梯度與高度的線性相關(guān)系數(shù)(圖3)。在87個(gè)點(diǎn)位上，線性相關(guān)系數(shù)的絕對(duì)值大于0.85，重力異常垂直梯度在垂向上與高度呈現(xiàn)出極強(qiáng)的線性關(guān)系。

圖3 模擬點(diǎn)位上Δgr與高度的線性相關(guān)系數(shù)Fig.3 Linear correlation coefficients between Δgr and height at each test point

在垂向上，Δgr的變化整體較為緩慢。表1為各個(gè)點(diǎn)位上高度7 500 m與5 000 m、7 500 m與2 500 m以及5 000 m與2 500 m的Δgr差值統(tǒng)計(jì)。可以看出，三者Δgr差異的絕對(duì)值平均值分別為3.3 E、8.1 E、4.8 E。按5 000 m延拓高差計(jì)算，由梯度近似造成的相應(yīng)的誤差絕對(duì)值平均值分別為1.7 mGal、4.0 mGal以及2.4 mGal，對(duì)當(dāng)前的航空重力觀測(cè)精度以及向下延拓的精度而言，這個(gè)誤差在可接受范圍之內(nèi)。

表1 不同高度面Δgr模型值的差異Tab.1 Differences of Δgr on different heights

為了進(jìn)一步說明前文中重力異常垂直梯度呈強(qiáng)線性分布的結(jié)論，采用線性外推計(jì)算2 500 m高度處的Δgr：

圖4 Δgr線性外推示意圖Fig.4 Demonstration of Δgr linear extrapolation

(3)

1.2 重力異常向上延拓與垂直梯度

球面上某一點(diǎn)重力異常的垂直梯度可以用該高度面上的重力異常表示：

(4)

利用重力異常進(jìn)行重力異常向上延拓的積分公式為：

(5)

式中，P′點(diǎn)為與式(4)中P點(diǎn)同一向徑且高出P點(diǎn)H的點(diǎn)。式(5)與式(4)在數(shù)學(xué)上是等價(jià)的，證明如下。

在離散形式下通常使用式(6)進(jìn)行重力異常向上延拓，以避免中心區(qū)域?qū)е路e分結(jié)果發(fā)散[14]:

(6)

在圖4中，由中值定理有：

(7)

式中，ΔgP2H可由重力異常向上延拓穩(wěn)定得到。

本文利用重力異常垂直梯度進(jìn)行延拓的過程分為3步：

1)根據(jù)向下延拓高度H，將航空重力異常向上延拓H，利用式(7)得到航空重力高度面上方某點(diǎn)的垂直梯度；

2)近似計(jì)算航空重力高度面的垂直梯度。同樣采用向上延拓后差分的方式，首先將重力異常向上延拓1/5H(1/5H為本文差分近似的經(jīng)驗(yàn)值)，然后差分得到垂直梯度近似值；

3)利用式(3)外推航空重力高度面下方某點(diǎn)的垂直梯度，代入式(2)進(jìn)行向下延拓。

2 臺(tái)灣地區(qū)航空重力異常向下延拓

本文在臺(tái)灣地區(qū)進(jìn)行航空重力向下延拓實(shí)驗(yàn)。為說明本文方法的有效性與實(shí)用性，首先采用EGM2008模型模擬觀測(cè)值進(jìn)行模擬實(shí)驗(yàn)[15]，然后對(duì)實(shí)測(cè)航空重力數(shù)據(jù)進(jìn)行處理。臺(tái)灣地區(qū)有豐富的航空重力觀測(cè)與地面重力測(cè)量數(shù)據(jù)，本文采用5 156 m高度的航空重力數(shù)據(jù)進(jìn)行向下延拓，以及3 360個(gè)地面重力觀測(cè)點(diǎn)進(jìn)行檢核[8]。圖5為臺(tái)灣地區(qū)實(shí)測(cè)航空重力格網(wǎng)值，圖6為該地區(qū)地面重力值分布及相應(yīng)重力異常。

圖5 臺(tái)灣地區(qū)航空重力格網(wǎng)值Fig.5 Airborne gravity anomaly grid values in Taiwan area

圖6 地面重力觀測(cè)點(diǎn)位分布及觀測(cè)值大小Fig.6 Ground gravity points and their gravity anomaly

考慮到臺(tái)灣地區(qū)地形復(fù)雜，使用GRAVSOFT軟件對(duì)航空重力觀測(cè)值及地面重力值進(jìn)行以15′ 分辨率為平滑半徑的殘差地形改正[16-17]。考慮到邊緣效應(yīng)，將航空重力在22°～25°N、120°～122°E區(qū)域內(nèi)進(jìn)行格網(wǎng)化，然后向下延拓，格網(wǎng)分辨率為2′。航空重力異常最大值為253.72 mGal，最小值為-159 mGal。在該區(qū)域內(nèi)有地面重力觀測(cè)值3 360個(gè)。地面重力異常最低點(diǎn)高程為0，最高點(diǎn)高程為3 942.61 m；重力異常最小值為-52.54 mGal，最大值為271.95 mGal，標(biāo)準(zhǔn)差為78.50 mGal。在向下延拓實(shí)驗(yàn)中，將航空重力異常分別向下延拓至各地面重力點(diǎn)，再將延拓值與各點(diǎn)重力觀測(cè)值進(jìn)行比較。

為了說明本文方法的有效性，首先使用同樣位置分布的EGM2008模型的模擬值進(jìn)行延拓，向下延拓精度見表2。

表2 EGM2008模擬數(shù)據(jù)向下延拓精度統(tǒng)計(jì)Tab.2 Statistics of downward continuation results of EGM2008 simulated gravity anomaly

表3 航空重力向下延拓精度統(tǒng)計(jì)Tab.3 Statistics of downward continuation results of airborne gravity observations

3 結(jié) 語

本文聯(lián)合拉格朗日中值定理以及重力異常垂直梯度徑向分布的特點(diǎn)，改進(jìn)了利用重力異常垂直梯度進(jìn)行航空重力異常向下延拓的方法，利用向上延拓后進(jìn)行差分的方法獲取重力異常垂直梯度，然后外推進(jìn)行向下延拓。主要結(jié)論如下：

1)與常用的泰勒級(jí)數(shù)一階展開方式相比，使用拉格朗日中值定理進(jìn)行垂直梯度外推然后進(jìn)行向下延拓具有更高的理論嚴(yán)密性。

2)重力異常垂直梯度在臺(tái)灣地區(qū)在海平面到海平面以上10 000 m的空間范圍內(nèi)，在垂向分布上與高度具有強(qiáng)線性關(guān)系。在本文進(jìn)行模擬實(shí)驗(yàn)隨機(jī)選取的100個(gè)點(diǎn)位中，87個(gè)點(diǎn)位的重力異常垂直梯度與高度的線性相關(guān)系數(shù)絕對(duì)值超過了0.85。

3)在使用EGM2008模擬實(shí)驗(yàn)中，直接模擬航空重力高度面上的垂直梯度，利用外推得到航空重力高度面以下的垂直梯度，然后進(jìn)行向下延拓，精度達(dá)到2.3 mGal，較僅利用航空重力高度面垂直梯度或向上延拓差分梯度進(jìn)行向下延拓精度有較大提升。

4)在差分計(jì)算航空重力高度面垂直梯度中，航空重力高度面上垂直梯度計(jì)算的不準(zhǔn)確會(huì)導(dǎo)致利用外推梯度延拓的結(jié)果精度反而下降，這是理論嚴(yán)密與數(shù)值精度之間的矛盾。將航空重力向上延拓相應(yīng)點(diǎn)位處的向下延拓高差，將向上延拓前后重力異常與高度進(jìn)行差分得到差分梯度值，僅利用該差分梯度進(jìn)行向下延拓，精度達(dá)到10.1 mGal，該結(jié)果亦優(yōu)于傳統(tǒng)延拓結(jié)果[8,19]。這也說明了使用垂直梯度進(jìn)行向下延拓的優(yōu)越性。

可以預(yù)期，在能得到航空重力高度面上較高精度垂直梯度的情況下，使用本文方法將進(jìn)一步提升向下延拓的精度。下一步研究將集中在優(yōu)化利用重力異常計(jì)算重力異常垂直梯度的方法上。