航空重力數據向下延拓的逆泊松積分算法

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

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

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

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

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

1 基本原理

1.1 逆泊松積分

假設函數在球外滿足拉普拉斯方程ΔV＝0，則調和函數V在球面上和外部空間的關系可表示為：

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

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

傅里葉變換具有時域卷積性質，即

其中，F表示二維傅里葉變換算子。根據時域卷積的定義，式（2）可表達為：

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

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

向下延拓可表示為：

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

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

式中，S（u，v）是一個濾波器。本文選取高斯濾波器和維納濾波器，通過實驗選擇最優濾波參數，并對比兩種結果。

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

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

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

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

無論選用何種濾波器，a都影響向下延拓的結果。a值的選擇必須經過嚴格的實驗選出最優值，通過調節參數a，可在噪聲與信號之間取得折衷。如沒有外部數據檢核，參數的選擇可根據目標的分辨率和精度，參照EGM08模型確定。

1.2 迭代法

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

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

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

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

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

2 數值計算和結果分析

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

在計算中應用了100%零填充技術以減弱循環卷積誤差和邊緣效應。圖2表明了100%零填充的原理，除去重力格網外其余全部設置為0。

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

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

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

選用FFT 算法向下延拓時，首先要選擇濾波參數，假設a＝5，10，15，…，根據計算得到的結果選擇最優值。表2為濾波參數與標準偏差的關系統計信息。

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

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

表3 最優濾波參數下重力異常向下延拓誤差統計Tab.3 Statistics of continued errors of gravity anomaly in optimal filter parameters

采用迭代法進行向下延拓，當迭代次數等于50時，計算值和實測值的差值小于閾值10－2μGal。表4給出了迭代法向下延拓的結果。

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

比較這兩種結果，快速傅里葉變換算法的標準偏差為4.15 mGal，優于迭代法。圖3 給出向下延拓結果與地面值差值的空間分布。

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

3 結 語

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

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

［1］王興濤，石磐，朱非洲.航空重力測量數據向下延拓的正則化算法及其譜分解［J］.測繪學報，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，Forsberg R.Combining Airborne and Ground Gravity Using Collocation［A］／／Geodesy on the Move［M］.Springer Berlin Heidelberg，1998

［3］孫中苗.利用空中平均重力異常確定區域大地水準面［J］.鄭州測繪學院學報，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］.測繪學報，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］.地球物理學報，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］蔣濤，李建成，黨亞民，等.基于矩諧分析的區域重力場建模［J］.中國科學：地球科學，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］孫中苗，夏哲仁，王興濤.利用航空重力確定局部大地水準面的精度分析［J］.武漢大學學報：信息科學版，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］孫文，吳曉平，王慶賓，等.航空重力數據向下延拓的波數域迭代Tikhonov正則化方法［J］.測繪學報，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，Forsberg R.The Use of FFT Techniques in Physical Geodesy［J］.Geophysical Journal International，1990，100（3）：485－514

［10］劉曉剛，李姍姍，吳星.衛星重力梯度數據的向下延拓［J］.大地測量與地球動力學，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］.地球物理學報，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）