航空重力數(shù)據(jù)向下延拓的逆泊松積分算法

盧雪盈 柳林濤 梁星輝 劉金釗

1 中國科學(xué)院測量與地球物理研究所，武漢市徐東大街340號，430077

2 中國科學(xué)院大學(xué)，北京市玉泉路甲19號，100049

航空重力數(shù)據(jù)向下延拓基本的解題思路是求解Poisson方程，即逆Poisson積分方法。向下延拓是一個不穩(wěn)定的過程，高頻噪聲經(jīng)向下延拓后被嚴(yán)重放大甚至淹沒重力信號本身。為提高向下延拓的精度，許多學(xué)者提出了不同的延拓方法。Tikhonov在1977年給出了求解不適定方程的Tikhonov 正則化法，此后該方法得到廣泛應(yīng)用［1－9］。這些方法在一定程度上改善了向下延拓的結(jié)果，但是也都存在一些問題。例如，最小二乘配置法需要已知可靠的先驗(yàn)協(xié)方差模型，而重力異常的協(xié)方差模型在實(shí)際中很難精確求得；正則化方法需要選取合適的正則化參數(shù)，若參數(shù)選擇不當(dāng)，則會引起很大的誤差；僅使用多尺度邊緣約束法，對向下延拓結(jié)果的改善非常有限，它只是對現(xiàn)有方法的補(bǔ)充。

本文基于逆泊松積分進(jìn)行求解，采用快速傅里葉變換和迭代法來計算重力異常。快速傅里葉變換法在計算重力場卷積積分時表現(xiàn)出很大的優(yōu)勢，并且計算速度快，原理相對較簡單，特別適合處理規(guī)則格網(wǎng)數(shù)據(jù)，易于編程。

1 基本原理

1.1 逆泊松積分

假設(shè)函數(shù)在球外滿足拉普拉斯方程ΔV＝0，則調(diào)和函數(shù)V在球面上和外部空間的關(guān)系可表示為：

其中，Vρ和VR分別為球外和球面的調(diào)和函數(shù)，ρ為球外一點(diǎn)到球心的距離，R為球平均半徑，r為球外點(diǎn)到積分元dω之間的距離，此即泊松積分。假設(shè)T為地球重力場擾動位，那么重力場的延拓在數(shù)學(xué)上可表示為大地測量第一邊值問題。航空重力異常向下延拓即是求解泊松積分方程，稱為逆泊松積分。

采用球面的平面近似，式（1）變換為：

傅里葉變換具有時域卷積性質(zhì)，即

其中，F(xiàn)表示二維傅里葉變換算子。根據(jù)時域卷積的定義，式（2）可表達(dá)為：

因此，式（2）在頻率域表示為：

其中，G（u，v）｜z＝h和G（u，v）｜z＝0分別為h處和海平面處重力值的傅里葉變換，此為向上延拓。圖1（a）給出了向上延拓的頻譜特征，隨著頻率的增加而快速下降。

向下延拓可表示為：

圖1（b）為向下延拓的頻譜特征，e2πhf這部分會放大高頻噪聲，很小的誤差即會引起很大的反應(yīng)。所以向下延拓是不穩(wěn)定的，需要在算法中引入低通濾波。設(shè)計的濾波器不僅要具有良好的窄帶低通特性，還要有尖銳的高頻衰減性能：

圖1 延拓的一維傳遞函數(shù)Fig.1 The 1Dtransfer function of continuation

式中，S（u，v）是一個濾波器。本文選取高斯濾波器和維納濾波器，通過實(shí)驗(yàn)選擇最優(yōu)濾波參數(shù)，并對比兩種結(jié)果。

高斯函數(shù)具有旋轉(zhuǎn)對稱性，即濾波器在各個方向上的平滑程度是相同的，其可分離性使得較大尺寸的高斯濾波器可以有效實(shí)現(xiàn)。高斯濾波的傳遞函數(shù)如下：

其中，a為濾波參數(shù)，a越大，平滑程度就越好。采用高斯濾波向下延拓表示為：

維納濾波器（Wiener filter）是一種以最小平方為最優(yōu)準(zhǔn)則的線性濾波器。在一定的約束條件下，其輸出與期望輸出的差的平方達(dá)到最小。在頻域，維納濾波器可表達(dá)為［11］：

式中，a＞0為濾波參數(shù)，同樣決定濾波器的平滑程度。利用維納濾波器平滑重力異常向下延拓表示為：

無論選用何種濾波器，a都影響向下延拓的結(jié)果。a值的選擇必須經(jīng)過嚴(yán)格的實(shí)驗(yàn)選出最優(yōu)值，通過調(diào)節(jié)參數(shù)a，可在噪聲與信號之間取得折衷。如沒有外部數(shù)據(jù)檢核，參數(shù)的選擇可根據(jù)目標(biāo)的分辨率和精度，參照EGM08模型確定。

1.2 迭代法

在實(shí)際計算中，迭代法是采用逐步趨近的方法進(jìn)行計算的。已知航線高度h處的重力異常，且大地水準(zhǔn)面外沒有質(zhì)量，迭代法向下延拓的步驟如下［12］：

1）將航線高度h處的重力異常Δg（x）垂直放置在h＝0，即地面的對應(yīng)點(diǎn)上，作為初值Δg（0）；

2）將Δg（0）利用公式（6）向上延拓到h處，得到G（1）；

3）根據(jù)h處實(shí)測值和G（1）的差異，對Δg（0）進(jìn)行校正，作為h＝0面的新重力異常值，然后以校正后的值再次進(jìn)行向上延拓計算，Δg（n＋1）＝Δg（n）＋g（x）－G（n）；

4）重復(fù)步驟2）、3），直到h處的實(shí)測值與G（n）的差值小于設(shè)定的閾值。

2 數(shù)值計算和結(jié)果分析

為評定航空重力數(shù)據(jù)的精度以及逆泊松積分法向下延拓的精度，采用某地區(qū)2 000m 高度的航空重力異常和地面實(shí)測重力異常，格網(wǎng)間距為100m×100m。航空重力測量的平均飛行速度為248km／h，北－南方向共有111 條每200 m 間隔的方向線（traverse line），東－西方向共10條每2 000m 間隔的控制線（control line）。將空中重力異常向下延拓到地面，與地面已知重力異常作比較，向下延拓采用FFT 方法和迭代法，積分半徑為30′。

在計算中應(yīng)用了100%零填充技術(shù)以減弱循環(huán)卷積誤差和邊緣效應(yīng)。圖2表明了100%零填充的原理，除去重力格網(wǎng)外其余全部設(shè)置為0。

圖2 延拓時采用的100%零填充原理Fig.2 The theory of zero－padding（100%）in continuation

評定航空重力數(shù)據(jù)的精度時，將地面已知重力異常按照泊松積分向上延拓到航線高度處，與空中重力異常進(jìn)行比較，其差值的標(biāo)準(zhǔn)差為3.58 mGal，具體結(jié)果如表1。

表1 空中重力異常的精度Tab.1 The accuracy of airborne gravity anomaly

選用FFT 算法向下延拓時，首先要選擇濾波參數(shù)，假設(shè)a＝5，10，15，…，根據(jù)計算得到的結(jié)果選擇最優(yōu)值。表2為濾波參數(shù)與標(biāo)準(zhǔn)偏差的關(guān)系統(tǒng)計信息。

由表2知，高斯濾波和維納濾波的參數(shù)分別為10和15時，延拓結(jié)果較好，于是在這兩個參數(shù)附近取值進(jìn)行計算，最后得出高斯濾波和維納濾波的最優(yōu)參數(shù)分別為10和13（表3）。

表2 重力異常向下延拓的誤差統(tǒng)計Tab.2 Statistics of continued errors of gravity anomaly

表3 最優(yōu)濾波參數(shù)下重力異常向下延拓誤差統(tǒng)計Tab.3 Statistics of continued errors of gravity anomaly in optimal filter parameters

采用迭代法進(jìn)行向下延拓，當(dāng)?shù)螖?shù)等于50時，計算值和實(shí)測值的差值小于閾值10－2μGal。表4給出了迭代法向下延拓的結(jié)果。

表4 迭代法重力異常向下延拓統(tǒng)計Tab.4 Statistics of continued errors of gravity anomaly in iterative method

比較這兩種結(jié)果，快速傅里葉變換算法的標(biāo)準(zhǔn)偏差為4.15 mGal，優(yōu)于迭代法。圖3 給出向下延拓結(jié)果與地面值差值的空間分布。

圖3 兩種方法向下延拓后延拓值和地面值差異的空間分布Fig.3 The spatial distribution of differences between continued value and ground value after continuation with these two methods

3 結(jié) 語

采用迭代法計算時，相鄰的兩次迭代值具有密切關(guān)系，前一次迭代存在的誤差，會影響下一次迭代的結(jié)果，反復(fù)迭代后相當(dāng)于誤差的累加。由于存在邊界效應(yīng)，每次迭代都損失一部分?jǐn)?shù)據(jù)，這使得測區(qū)范圍不大時，在實(shí)際應(yīng)用中難以達(dá)到要求。FFT 算法運(yùn)算速度快，原理簡單，但是由于采用的是平面近似公式，當(dāng)延拓高度較大時，精度會降低，而且忽略了遠(yuǎn)區(qū)貢獻(xiàn)的影響引起的誤差；解卷積過程有自身的缺點(diǎn)，e2πhf對數(shù)據(jù)中的高頻噪聲起放大作用，很小的一個誤差即會對最后結(jié)果產(chǎn)生巨大的影響，因此在應(yīng)用該方法時，應(yīng)引入低通濾波。

本文利用實(shí)測航空重力數(shù)據(jù)，通過對FFT 算法和迭代法延拓結(jié)果進(jìn)行對比和分析，評價以上兩種方法的精度和可靠性。在航空重力異常精度為3.58mGal的情況下，利用FFT 算法引入低通濾波后，計算得到高斯濾波和維納濾波最優(yōu)參數(shù)分別為10 和13，相應(yīng)的延拓精度為4.31 mGal和4.15 mGal；利用迭代法計算的延拓精度為4.97mGal。兩種延拓結(jié)果都表明，向下延拓是一個不穩(wěn)定的過程。在不考慮地形影響的前提下，F(xiàn)FT 算法結(jié)果較好；引入低通濾波后，在都采用最優(yōu)濾波參數(shù)的情況下，維納濾波器的性能要優(yōu)于高斯濾波器。

［1］王興濤，石磐，朱非洲.航空重力測量數(shù)據(jù)向下延拓的正則化算法及其譜分解［J］.測繪學(xué)報，2004，33（1）：33－38（Wang Xingtao，Shi Pan，Zhu Feizhou.Regularization Methods and Spectral Decomposition for the Downward Continuation of Airborne Gravity Data［J］.Acta Geodaetica et Cartographica Sinica，2004，33（1）：33－38）

［2］Tscherning C C，Rubek F，F(xiàn)orsberg R.Combining Airborne and Ground Gravity Using Collocation［A］／／Geodesy on the Move［M］.Springer Berlin Heidelberg，1998

［3］孫中苗.利用空中平均重力異常確定區(qū)域大地水準(zhǔn)面［J］.鄭州測繪學(xué)院學(xué)報，1999，16（4）：241－243（Sun Zhongmiao.The Determination of Regional Geoid with Aerial Mean－Gravity Anomalies［J］.Journal of the Zhengzhou Institute of Surveying and Mapping，1999，16（4）：241－243）

［4］石磐，王興濤.利用航空重力測量和DEM 確定地面重力場［J］.測繪學(xué)報，1997，26（2）：117－121（Shi Pan，Wang Xingtao.Determination of the Terrain Surface Gravity Field Using Airborne Gravimetry and DEM［J］.Acta Geodaetica et Cartographica Sinica，1997，26（2）：117－121）

［5］寧津生，汪海洪，羅志才.基于多尺度邊緣約束的重力場信號的向下延拓［J］.地球物理學(xué)報，2005，48（1）：63－68（Ning Jinsheng，Wang Haihong，Luo Zhicai.Downward Continuation of Gravity Signals Based on the Multiscale Edge Constraint［J］.Chinese Journal of Geophysics，2005，48（1）：63－68）

［6］蔣濤，李建成，黨亞民，等.基于矩諧分析的區(qū)域重力場建模［J］.中國科學(xué)：地球科學(xué)，2014，44（1）：82－89（Jiang Tao，Li Jiancheng，Dang Yamin，et al.Regional Gravity Field Modeling Based on Rectangular Harmonic Analysis［J］.Science China：Earth Sciences，2014，44（1）：82－89）

［7］孫中苗，夏哲仁，王興濤.利用航空重力確定局部大地水準(zhǔn)面的精度分析［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2007，32（8）：692－695（Sun Zhongmiao，Xia Zheren，Wang Xingtao.Determining Accuracy of Local Geoid Determined from Airborne Gravity Data［J］.Geomatics and Information Science of Wuhan University，2007，32（8）：692－695）

［8］孫文，吳曉平，王慶賓，等.航空重力數(shù)據(jù)向下延拓的波數(shù)域迭代Tikhonov正則化方法［J］.測繪學(xué)報，2014，43（6）：566－574（Sun Wen，Wu Xiaoping，Wang Qingbin，et al.Wave Number Domain Iterative Tikhonov Regularization Method for Downward Continuation of Airborne Gravity Data［J］.Acta Geodaetica et Cartographica Sinica，2014，43（6）：566－574）

［9］Schwarz K P，Sideris M G，F(xiàn)orsberg R.The Use of FFT Techniques in Physical Geodesy［J］.Geophysical Journal International，1990，100（3）：485－514

［10］劉曉剛，李姍姍，吳星.衛(wèi)星重力梯度數(shù)據(jù)的向下延拓［J］.大地測量與地球動力學(xué)，2011，31（1）：132－137（Liu Xiaogang，Li Shanshan，Wu Xing.Downward Continuation of Satellite Gravity Gradient Data［J］.Journal of Geodesy and Geodynamics，2011，31（1）：132－137）

［11］Forsberg R，Skourup H.Arctic Ocean Gravity，Geoid and Sea－Ice Freeboard Heights from ICESat and GRACE［J］.Geophysical Research Letters，2005，32（21）

［12］徐世浙.迭代法與FFT 法位場向下延拓效果的比較［J］.地球物理學(xué)報，2007，50（1）：285－289（Xu Shizhe.A Comparison of Effects between the Iterative Method and FFT for Downward Continuation of Potential Fields［J］.Chinese Journal of Geophysics，2007，50（1）：285－289）