鄭 崴，陳學鋒

（1.中國地質大學 （北京）地球物理與信息技術學院，北京 100083；2.河南高等機電?？茖W校，河南 新鄉(xiāng) 453000；3.河南科技學院信息工程學院，河南 新鄉(xiāng) 453003）

1 緒言

重力勘探是地球物理勘探的基本方法之一。它是在地球表面上利用重力儀，測量由于地下物質密度分布不均勻引起的重力變化（或重力異常），這種物質密度的分布與地質構造及礦產分布有密切的關系。因此，通過研究重力異常，可以了解和推斷地球的結構、地殼的構造，以及勘探礦產資源等等。最早用于重力勘探的設備是測量重力梯度的扭稱，但是由于設備笨重，效率低下，在20世紀50年代，被更輕便，效率較高的重力儀所取代［1］?？梢哉f重力勘探是由重力梯度測量開始的。20世紀70年代以后，出于軍事方面應用的需求，重力梯度儀獲得了飛速的發(fā)展，美國、法國等國家先后研制出了不同原理、不同結構的重力梯度儀。而現(xiàn)代電子技術、計算機技術、超導量子干涉技術、低溫微波空腔諧振技術、超導負彈簧技術的發(fā)展和應用，重力梯度儀的靈敏度和穩(wěn)定性取得了突破性的進展。由于重力梯度相對于重力異常具有更高的分辨率，隨著重力梯度儀的發(fā)展，它在地質資源勘查、地質調查等領域發(fā)揮出重要的作用。

相對于地面測量技術來講，航空重力梯度測量具有快速、經濟和高精度的特點，它幾乎可以到達任何測量區(qū)域。除了應用于精細地質勘探外，還應用于國防、軍事領域。因此，航空重力梯度測量技術是重力勘探中最新、最先進的技術。本文首先介紹航空重力梯度測量的基本概念與發(fā)展現(xiàn)狀，然后討論標準傅里葉變換在航空重力梯度數據預處理中的應用。

2 重力和重力梯度

在地球表面任一點的重力是由地球引力與地球自轉產生的慣性離心力的合力，它與在地球表面的位置和物體的質量有關。三度體在u＝（x，y，z）處的重力場可以用重力位的標量函數

式中，G為萬有引力常數，為u'＝（x'，y'，z'）物質的密度。

在笛卡爾坐標系中，可以從重力位推導出重力和重力梯度的表達式。重力是重力位的一階導數，即

重力儀測量的是垂直方向分量gz，即

而重力梯度是重力位的二階導數，即重力加速度矢量的三個分量分別沿直角坐標系三個方向的一階導數。重力梯度是一個張量，可以表示為式（4）中的形式。圖1位重力場矢量及其梯度示意圖。

在無源的空間內，引力的散度和旋度均為0，所以重力梯度是一個對稱的張量，有TXY＝TYX，TXZ＝TZX，TYZ＝TZY，且有TXX＋TYY＋TZZ＝0，因此，重力梯度張量的9個分量中，只有5個是獨立的分量。

在球坐標系中，設有一質量為M的均質球體，則有

圖1 重力場矢量及其梯度示意圖

式中，r為測量點到目標體中心的距離，由式（5）可見，物體的密度與距離r成反比，重力與r2成反比而重力梯度與距離的3次方r3成反比，可見，重力梯度比重力衰減的快，利用重力梯度可以更好的估計物體邊緣的位置，所以說重力梯度相比于重力來講具有更高的分辨率。而重力梯度異常值反映的是由于地下密度分布不均勻的重力異常在空間的變化率，重力梯度張量的各分量，可以反映出目標地質體不同的地質特征，如圖2（a）是針對已知地質目標完成的一次航空重力梯度測量所獲取的結果，已知目標地質體為鹽類的蓋層，表現(xiàn)為TZZ圖中圓圈部分所顯示的正值異常。而在圖2（b）中的各個張量分量則可以明顯的分辨出目標體的地質特征，其中TXX和TXY可以顯示出目標體的南北向的邊緣位置，從圖中能夠明顯的分辨出目標體的南北向的位置輪廓；同樣的，由TYY和TYZ圖可以分辨出目標體的東西向的邊緣位置。而由TXY圖則顯示出蓋層具有典型的2高點和2低點的特點4極異常。

圖2 已知鹽類蓋層航空重力梯度測量張量分量特征

3 航空重力測量的基本原理和航空重力梯度測量系統(tǒng)

根據牛頓第二定律，作用于單位質點上的總加速度矢量fi（稱為比力），與載體運動加速度矢量ri和引力加速度矢量Gi的關系可以表示為

因此，要獲得引力加速度Gi，必須獲得載體運動加速度ri，再從測得的比力fi中將其減去。在地面靜止狀態(tài)下，重力儀的運動加速度ri＝0，此時重力儀測得的就是重力場矢量。而在航空重力梯度測量中，由于運動加速度不斷的變化且不為零，就必須在測得的比力fi中分離出運動加速度ri，而分離方法是通過對共用基線的兩個加速度計的輸出求差，以消除載體運動的影響。如果共用基線是旋轉穩(wěn)定的，由差值讀數可以獲得重力梯度分量，這也是重力梯度測量的基本原理。除了要消除載體運動加速度帶來的影響，還要保證在測量過程當中，加速度計保持在平穩(wěn)、穩(wěn)定的工作狀態(tài)，這一點主要是通過硬件來實現(xiàn)，如可以將加速度計安裝在INS（Inertial Navigation System）慣導平臺等。

這里主要介紹兩種已經投入商用的航空重力梯度測量系統(tǒng)。

3.1 基于旋轉加速度計的重力梯度系統(tǒng)

旋轉加速度計重力梯度系統(tǒng)是市面上唯一投入商用的航空重力梯度系統(tǒng)，又分為全張量測量系統(tǒng)和部分張量測量系統(tǒng)。

全張量重力梯度測量系統(tǒng)是由美國Bell Aerospace（現(xiàn)已并入Lockheed Martin公司）自主研究制造的Air-FTG系統(tǒng)。它的核心器件都是圖3所示的旋轉加速度儀（GGI，Gravity Gradient Instrument），它是將兩對相互正交的加速度計安裝在一個圓底盤的邊緣上。GGI旋轉工作時，每對加速度計感應的圓盤的運動加速度大小相等方向相反，可以互相抵消。通過計算每對加速度計間的距離和底盤的旋轉頻率，可以測得一個平面上的兩個梯度分量。一個全張量梯度測量系統(tǒng)包括三個GGI，這三個GGI的軸線互相垂直，每個GGI的軸線沿鉛垂線以同樣的角度相交，從上向下看，三條軸線的投影角度為120度，如圖4所示。為保證測量精度，GGI必須以恒定的速度繞垂直軸方向旋轉。同時，為了隔離載體運動過程對系統(tǒng)的影響，可以將旋轉加速度儀安裝在由陀螺儀穩(wěn)定的慣性平臺上。另外，系統(tǒng)配備獨立的坐標系，通過將三個GGI在外部坐標上的梯度分量進行線性組合，可以計算出重力梯度張量的5個獨立分量，進而得到重力梯度張量的全部分量值。

圖3 加速度計安裝及轉盤運動示意圖

圖4 Air-FTG的三軸GGI示意圖

部分重力梯度張量測量系統(tǒng)的代表是由澳大利亞的BHB Billiton公司出資研制的FALCON系統(tǒng)，此系統(tǒng)于1997年研制成功，1999年投入商業(yè)航空物理勘探。與全張量系統(tǒng)相比，F(xiàn)ALCON只配備一個 GGI，它可以直接測量出TXY，（TXX－TYY）/2兩個梯度分量，再通過面積性測量和數學變換得到重力水平分量gx和梯度分量TXX，根據實際需要，它還可以由TXY計算出其他的張量分量，也就是說FALCON系統(tǒng)可以通過數學計算的方法得到全部張量分量。

3.2 超導航空重力梯度系統(tǒng)

超導航空重力梯度系統(tǒng)（SGG）的研制源于20世紀70～80年代，近十年來，超導航空重力系統(tǒng)作為下一代的主力航空重力勘探系統(tǒng)，已經成為目前研究的重點和熱點。超導航空重力梯度系統(tǒng)的設計指標是具有更高的靈敏度和更精確的測量結果，設計噪聲指標優(yōu)于1E。

英國的ARKeX公司正在改造用于衛(wèi)星上超導重力梯度系統(tǒng)，由于衛(wèi)星的飛行比較平穩(wěn)，而飛機的加速度很大，ARKeX改造的重點就在克服飛機的飛行干擾上。ARKeX公司的最終目標是生產經營世界上首臺超導勘探重力梯度儀EGG（Exploration Gravity Gradiometer）。EGG的設計可配備六個加速度計模塊，用來測量全部的梯度張量分量，但在實驗過程中制造出的第一臺EGG只用了兩個加速度計模塊，只能測量垂直方向的張量元素TZZ。

4 航空重力梯度測量的數據預處理

全張量重力梯度測量系統(tǒng)是目前應用較廣泛的航空重力梯度測量系統(tǒng)，以它的數據處理過程來介紹航空重力梯度的數據預處理方法。FTG的數據處理過程主要有兩個階段：數據預處理階段和后任務處理階段。第一個階段是實現(xiàn)對加速度計的輸出值因外界測量環(huán)境的變化而進行的補償。加速度計是在飛機高速運動的狀態(tài)下完成測量，因此，在這個階段，需要對量測量進行一些必要的校正，這些校正主要針對以下幾種影響：①每一個AGG的底部磁盤高速旋轉產生的向心力的影響；②飛機飛行方向變化和飛行姿態(tài)變化的影響；③飛機飛行過程中燃油消耗量變化的影響④飛行過程中加速度不斷變化的影響。

當這些校正完成后，還要對校正后的數據進行插值，以供第二個數據處理階段——后任務處理階段使用。在預處理的最后，需要將校正后的數據重新插值到0.5Hz，同時解算出每一個張量元素值，并將其換算到外部現(xiàn)實世界的直角坐標系統(tǒng)。而后任務處理階段的主要任務是降噪，即在存在大量隨機干擾的測量數據中獲得有用數據的信息，使得輸出的張量分量具有較為滿意的信噪比。這需要多種濾波技術的支持來實現(xiàn)。

4.1 傅里葉變換解算重力梯度張量分量

不難看出，在數據預處理階段，解算出重力梯度張量的各個分量并實現(xiàn)坐標軸的變換，在數據預處理階段是非常重要的，為了獲得符合解釋要求的數據，還要保證解算出的張量具有可觀的信噪比（S/N）。這種數據預處理的要求，可以由標準傅里葉變換的方法來實現(xiàn)。因為重力梯度張量實際是重力位的在各個坐標軸方向的二階導數，利用重力梯度的傅里葉變換就可以導出相應的重力位。重力位是重力測量中很重要的物理量，一旦確定了重力位，所有的重力分量都可以由它得到。并且，在進行傅里葉變換的同時，可以設計合適的加權函數來抑制噪聲的影響，達到提高信噪比的目的。

設在測量過程中飛機的飛行高度是不變的，在高度為Z的X-Y水平面內，重力位的傅里葉變換對為

對式（9）的傅里葉變換部分求積分，得

不難看出，重力的傅里葉變換和重力位的傅里葉變換存在著一定的關系，再經過反變換，可得其重力位

所以，在高度為Z水平面內，只要測出gx，gy，TXX，TYY，TXY五個量中的任意一個分量，就可以通過相應的變換關系得出重力位，此時，所有X，Y方向的重力的導數就都可以計算出來［11］。而在實際的數據處理過程中，設需要處理的FTG測量數據已經實現(xiàn)網格化，按照網格的水平方向和垂直方向對數據進行傅里葉變換，而垂直方向的傅里葉變換相對于水平方向的傅里葉變換更容易實現(xiàn)，此時式（7）的關系可表示為式（12）。

在進行傅里葉變換時，可以將整個待處理數據劃分為多個小的數據塊，每個小的數據塊又有多個網格組成。對這些由多個網格組成的數據進行積分運算，可推出相應的重力位，再求重力位的均值，可以提高輸出的信噪比，再進行微分變換，則可以推導出其他的張量分量，如TZZ。在這個計算過程中，為導出具有較高信噪比的重力分量的均值，可以通過將重力分量轉換到頻率域內，根據量測數據的頻譜結構，設計相應的加權函數濾除域內的高頻噪聲或者放大有用信號的功率，以此來達到提高信噪比的目的。

所以，在理論上來講，由重力位或其任意導數的傅里葉變換都可以導出梯度張量的各個分量，反之亦然。但在實際應用中，有些變換是無法直接實現(xiàn)的，比如說利用TXX求垂直方向的重力分量gz，理論上可以由式（14）直接計算得到。

但是，在式（14）中，由于TXX是沿X軸方向分量，它不包含垂直方向Z的信息，這就給數據的處理帶來兩個問題。第一個問題，式（14）中的分母kx≈0，即無法利用TXX的傅里葉變換求得重力分量gz。第二個問題，接近ky軸的地方，TXX的幅值衰減的很快，而噪聲信號的幅值卻很大，若在k域內直接進行傅里葉變換，相當于對噪聲信號進行放大，這必然導致信噪比的降低。

對于噪聲的影響，可以有多種方法來解決?？梢詫V波器的截止頻率設置在遠離ky軸的位置上，或者設計維納濾波器來減弱噪聲的影響。而對于第一個問題，常用的方法是，利用已測得的包含垂直方向的張量分量（如TXZ，TYZ）來導出gz，即進行傅里葉變換前對這些測得的張量分量進行線性組合。經過線性組合后，有用信號的頻譜會較為完整的保留下來，再去利用傅里葉變換去導出重力位，則可以計算出具有較高信噪比重力分量gz。但是這種方法的有兩方面的問題，首先，進行數據的線性組合需要完成大量的運算。其次，線性組合后會給降噪帶來一定的難度，如果有用的數據與噪聲數據組合在一起，數據變換后很難用濾波器去削弱這部分噪聲的影響。所以，利用線性組合來解決這個問題并不是最好的選擇。

而FTG系統(tǒng)則是用硬件設計的方式來解決這一問題，一個FTG系統(tǒng)包含三個GGI，每個GGI的軸線與垂直方向有相同的夾角120°。也就是說，每個GGI的輸出的張量分量當中都包含有關垂直方向的梯度元素TZZ的分量，因此，GGI輸出結果中包含垂直方向的信息。設三個GGI輸出的橫差曲率分別表示為C1，C2和C3，則梯度分量TZZ可以直接有公式（15）計算得出

4.2 坐標系的變換

實測數據是基于測量坐標系（X軸指向地理東，Y軸指向地理北，Z軸指向垂直方向），由于在傅里葉變換中要用到旋轉矩陣，這增加了數學計算的難度。為此，可以利用式（13）構建一個差分矩陣因子，實現(xiàn)由重力位到直角坐標系的梯度張量TXYZ的轉換

如果在球坐標系中，有一重力梯度張量TUVW，則定義一個矩陣R，將其由直角坐標系換算到球坐標系中，此時的變換公式表示為式（17）

舉個例子來講，設測量所得的梯度分量為Tuw，已知它是一個包含有橫差曲率的分量（類似前邊的C3），則它的傅里葉變換表示為式（18）

只有在坐標原點上，式（17）的分母為0。而丟失的數據部分為一平均值，在整個數據處理過程它是一個常量，因此不會影響后續(xù)的處理結果。

利用式（18）對TUVW的各個張量元素進行推導，可以導出相應的重力位，同時，傅里葉變換后，平均運算將有用信號平均分配到每個張量分量中，以此來實現(xiàn)達到提高有用信號功率，抑制噪聲的目的。

5 總結

地球重力梯度的研究一直是地球物理勘探領域最前沿、最熱點的研究課題，重力梯度實際上是重力異常在空間內各個方向上的變化率，可以更好的反映出探測目標的地質信息，相對于重力異常而言重力梯度異常具有更高的分辨率。同時，依托于航空平臺的重力梯度測量系統(tǒng)相對于傳統(tǒng)的地面測量儀器來說，又具備快速，高效，經濟的特點，具有廣闊的研究前景和價值。

由于其是在高速動態(tài)的狀態(tài)下完成實地的測量，航空重力梯度測量的數據處理步驟不同于傳統(tǒng)重力測量的處理方法，在數據預處理的過程中，利用傅里葉變換可以快速的對重力梯度數據進行預處理，以獲得相應的重力梯度張量分量和所需要的重力分量，同時也可以快速實現(xiàn)梯度張量在球坐標與直角坐標間的變換。

［1］曾華霖.重力場與重力勘探［M］.北京：地質出版社，2005.

［2］張金煥，李曉斌.航空重力梯度測量系統(tǒng)發(fā)展及現(xiàn)狀［J］.中國礦業(yè)，2011，20（zk）：188-201.

［3］張永明，張貴賓，盛君.航空重力梯度測量技術及應用［J］.工程地球物理學報，2006，3（5）：375-380.

［4］張永明，張貴賓，盛君.航空重力梯度測量的基本理論與應用［J］.地質裝備，2006，7（6）：23-27.

［5］張昌達.航空重力測量和航空重力梯度測量問題［J］.工程地球物理學報，2005，2（4）282-290.

［6］紀兵，劉敏，呂良，等.重力梯度儀在水下航行中的應用［J］.海洋測繪，2010，30（7）：23-25.

［7］孫中苗，翟振和.航空重力測量數據的小波閾值濾波［J］.武漢大學學報：信息科學版，2009，34（10）：1222-1225.

［8］Richard Lane.Airborne Gravity 2004：Abstract from the ACEG-PESA Airborne Gravity 2004workshop ［A］.Geoscience Australia，2004.

［9］Scott Hammond，Colm Murphy.Air-FTGTM：Bell Geospace's Airborne Gravity Gradiometer-A Description and Case Study.Preview，August 2003：24-26.

［10］王靜波，熊盛青，周錫華，等.航空重力測量系統(tǒng)研究進展［J］.物探與化探，2009，33（4）：368-373.

［11］張永明，張貴賓，盛君.航空重力梯度測量的基本理論及應用［J］.地質裝備，2006，7（6）：23-27.

［12］舒晴，周堅鑫，尹航，等.應用和研制中的航空重力梯度測量系統(tǒng)［C］.中國地球物理年會論文集，2011.

［13］Zhangxin Liu，Hai Lu.Research on Correlation Imaging of Gravity Gradient Based on FFT［C］.2011Xian International Conference on Fine Geological Exploration and Groundwater＆ Gas Hazards Control in Coal Mines，2011：203-209.

［14］Majid Beiki.Analytic Signals of Gravity Gradient Tensor and Their Application to Estimate Source Location ［J］.Geophysics，2010，75（6）：159-173.

［15］Gary Barnes，John Lumley.Processing Gravity Gradient Data［J］.Geophysics，2011，76（2）：133-147.

［16］Glennie，C.，and K.P.Schwarz.Airborne gravity by strapdown INS＝DPGS in a 100km by 100km area of the Rocky Mountains：Proceedings of the International Symposium on Kinematic Systems in Geodesy［J］.Geomatics and Navigation（KIS97），1997：619-624.

［17］Vassiliou，A.A.Stage II processing of airborne gravity gradiometer data using frequency domain techniques［C］.Proceedings of the 14th Annual Gravity Gradiometry Conference，1986：435-440.

［18］Dransfield，M.H.，Conforming Falcon gravity and the global gravity anomaly.Geophysical Prospecting，2010，58，469-483.

［19］Li Y.，and D.W.Oldenburg，Rapid construction of equivalent sources using wavelets［J］.Geophysics，2010，75（3），51-59.

［20］Saad，A.H.，Understanding gravity gradients—A tutorial［J］.The Leading Edge，2006（25），942-949.