孫中科，朱自強，魯光銀，曹書錦

（中南大學 地球科學與信息物理工程學院，湖南 長沙 410083）

0 前言

在地球物理勘探中，常規(guī)的重力勘探儀器只能測量到重力位的垂向一階導數(shù)，而重力梯度測量儀能夠測量重力位的二階導數(shù)，即重力梯度張量分量中獨立的五個分量［1］。由于重力張量擁有的各種優(yōu)異特征，國外對重力梯度張量特征的研究已經(jīng)進入到較深的層次［2～5］。Vasco［6］建立了一組長方體模型，在下底埋深不同的條件下，對重力張量的三個對角線元素（Vxx，Vyy，Vzz）的反演結果特征進行了分析；Zhdanov［7］基于三維正則化條件下的收斂技術對重力梯度的全部張量數(shù)據(jù)做了反演。這種方法的特點是在正演擬合過程中，采用的是同體積同質心的球作為基本剖分單元，一定程度上簡化了反演過程中的正演迭代計算。在對重力張量實測數(shù)據(jù)的去噪方面，Julio［8］首次提出將一維自適應引入小波濾波技術，這種方法在保留局部尖峰信號的同時，可以濾掉高頻噪聲；While［9］通過頻譜分析方法，得到了重力張量的二維功率及其法則，該法則可用于判斷所測得的重力梯度張量異常數(shù)據(jù)是否由地質原因引起；Kevin L.Mickus［10］通過對數(shù)據(jù)進行快速傅立葉變換（FFT），對垂向重力數(shù)據(jù)進行推導得到了全部的重力梯度張量數(shù)據(jù)。在此基礎上，在重力梯度張量解釋方面已有成功應用實例的還有功率譜反演方法［11］和歐拉褶積法［12］等。與國外相比，國內研究重力梯度張量的水平還處于比較落后的階段，且國內對重力梯度的研究主要集中在單獨計算某一個單獨的梯度張量分量上，在確定場源邊界上應用較多［13，14］，而在張量意義上的應用卻很少。魏偉［15］給出了一種能夠對斷裂的重力異常數(shù)據(jù)進行改進的梯度解釋方法；趙德軍［16］對水平重力梯度分量的邊值問題的級數(shù)解進行了推導；李姍姍［17］給出了基于重力矢量的Neumann邊值問題和Dirichlet邊值問題的解；劉尚余［18］在對衛(wèi)星重力梯度數(shù)據(jù)處理的過程中引入了快速傅里葉變換?？偟膩碚f，國內重力張量的研究還出于起步階段，但已經(jīng)展現(xiàn)出良好的發(fā)展勢頭。

粒子群優(yōu)化算法（簡稱PSO）是1995年Kenndey和Eberhart［19］提出的一種基于群體智能的優(yōu)化算法。粒子群優(yōu)化算法的特點是利用個體在解空間中的隨機速度來改變個體，其解群相對進化代數(shù)而言，標尺暗處更強的隨機性，其計算復雜度比其它算法低，搜索速度更快。吳招才等［20］將PSO算法應用于板狀體磁異常的反演，并與遺傳算法（GA）進行了比對；邱寧等［21］將混沌思想應用于PSO算法，利用混沌擾動的方法來解決PSO算法中的局部極值問題，并對磁法數(shù)據(jù)進行了反演，取得了一定的效果。

1 正演基礎

重力梯度張量被定義為重力位［22］在各個方向的二階導數(shù)。在笛卡爾坐標系條件下，假設重力位為V，則重力梯度張量可用一個3×3的矩陣來表示：

在無源空間內，引力的散度和旋度都為零，重力梯度張量矩陣為對稱矩陣，所以有：Vxx＋Vyy＋Vzz＝0，Vxy＝Vyx，Vxz＝Vzx，Vzy＝Vyz，可以看出，重力梯度的九個張量中只有五個是獨立的。重力梯度的國際標準單位為1／s2，但這個單位過于大，通常使用厄缶（E）作為重力梯度張量的量綱，1E＝1mGal／（10km）。

作者在本文中，將目標體劃分為許多小的長方體，計算每個長方體的張量后疊加，長方體模型如圖1所示。在圖1中，（ζ1，η1，ξ1）、（ζ2，η2，ξ2）分別為長方體角點坐標；（x，y，z）為測點坐標。

布格重力異常及重力各個張量分量計算公式如下［23、24］：

圖1 長方體正演模型Fig.1 The forward model of rectangular parallelepiped

2 粒子群優(yōu)化算法

粒子群優(yōu)化算法［21］（PSO）來源于對鳥群有序且不可預測的運動特點進行抽象化，將群體中的每一個個體都看作搜索空間（解空間）中某一個粒子，分別代表著反演問題的每一個可行解，每一個粒子具有速度和位置兩個基本屬性。每個粒子都依據(jù)自身及周圍粒子的迭代經(jīng)驗進行更新，即每一個粒子都通過評價自身最優(yōu)解和群體的最優(yōu)，來不斷對自己的速度和位置進行改正，粒子的優(yōu)劣程度（與反演模型的擬合程度）是通過對粒子所對應的適應度函數(shù)值確定的。粒子群算法求解過程不需要導數(shù)信息，只要給出的問題能夠求出適應度函數(shù)值，就能夠通過粒子自身的迭代進行求解，粒子群優(yōu)化算法基本的迭代公式如下：

其中 t是迭代次數(shù)；i是粒子在群體中的編號。

為迭代第t次時候粒子i的速度。

式（12）為粒子迭代第t次時候編號為i粒子所在位置；Pbest是粒子i的個體最優(yōu)位置；Gbest是解空間中全部粒子的最佳位置；r1和r2為［0，1］區(qū)間內均勻且獨立分布的任意數(shù)；c1和c2為迭代過程的學習因子，用來控制每次粒子更新時的步長；ω為慣性權重因子，用來控制粒子繼承上一代粒子特征的大小。式（9）是根據(jù)粒子自身位置、歷史最優(yōu)和全局最優(yōu)位置來計算獲得的最新速度。在迭代過程中，每一個粒子的速度都被限制在［－vmax，vmax］范圍內，用以控制粒子的極值不能超出解空間。然后群體根據(jù)式（10）飛行到新的位置。具體流程如下［25］：

（1）初始化參數(shù)，包括慣性因子ω、學習因子c1和c2。

（2）隨機產(chǎn)生M個粒子的種群，在允許范圍內隨機設置粒子的初始位置和速度。評價每一個粒子的適應值，由目標函數(shù)計算每一個粒子的適應度值并評價其優(yōu)劣。

（3）粒子在解空間搜索，對每一個粒子，其本身的最優(yōu)解記為Pbest。

（4）整個種群目前搜索到的最優(yōu)解稱為全局最優(yōu)解記為Gbest。

（5）每個粒子根據(jù)式（9）、式（10）來更新自身的速度和位置，并把速度限制在［－vmax，vmax］范圍內。

（6）檢查終止條件（達到設定迭代次數(shù)或最優(yōu)解滿足條件），如果滿足則終止迭代輸出結果，否則返回步驟（3）。

3 重力張量粒子群優(yōu)化反演

3.1 反演理論

地球物理反演問題可以表示為式（13）。

Gαβ＝Sαβ，σ （13）將Sαβ看作映射算子（即雅可比矩陣），σ為待測模型的密度參數(shù)，那么反演問題就是求已知觀測數(shù)據(jù)Gαβ下通過S－1αβ求反演模型σ的一種運算。在重力張量單分量反演中，分別對式（13）中的模型觀測數(shù)據(jù)和雅可比矩陣取相應的分量值即可。在包含全部張量的聯(lián)合反演中，我們選取其中五個獨立的分量（Vxx，Vxy，Vyy，Vxz，Vyz），令

式中 G為五個獨立分量的觀測值構成的矩陣；S為相應的幾何函數(shù)構成的矩陣。分別用式（14）、式（15）代替式（13）中的兩個矩陣，就得到重力全張量聯(lián)合反演的關系式。

3.2 目標函數(shù)

反演過程可以歸結為使如下最小二乘目標函數(shù)趨于極?。?/p>

式中 gi（i＝1，2，…，n）為n個離散反演模型的觀測值；gi（σ）（i＝1，2，…，n）為模型σ的理論正演場值在與gi點對應的n個離散采樣值；σ為反演迭代中的模型參數(shù)。

3.3 模型試算

為了驗證該算法的正確性，本文作者設計了如下目標體：將場源所在空間劃分成15×15×10個基本物性單元，每個基本單元在x、y方向上的長度都為40m，z方向上設為為50m，目標體模型為一個大小為200×200×200的正方體，其頂面埋深為200m，剩余密度設定為1×103kg／m3，其余點的密度為“0”。地表的測網(wǎng)共有15×15＝225個測點，沿x、y方向的點距都為40m。

反演參數(shù)設置如下：慣性權重因子ω取阻尼慣性權重，為（0，1）之間的某個數(shù)，在反演中可以變化；學習因子c1＝c2＝2；初始種群數(shù)設置為模型單元個數(shù)的二倍；初始解的取值范圍皆為（0，1）內的隨機數(shù)；反演中將粒子速度限制在（0，1）范圍內；迭代次數(shù)都設為500次。

對該目標體分別對各個張量單分量、布格重力異常和全張量進行反演，反演耗時約460s～490s（因初始解是隨機設置的，同一張量反演的耗時也會因為初始解設置的差別而略有不同），反演結果見圖2至圖5、及下頁圖6～圖9，并對比其效果。

從各個分量反演結果圖分析可以得到如下結論：

（1）Vxx、Vxy和Vxz分量反演結果都能夠大致體現(xiàn)出目標體的范圍，但總體上反演得到的密度值偏小，Vxy反演結果在形態(tài)上更接近初始模型。

（2）Vyy分量反演的效果較好，在密度值上較接近原始模型，但在形態(tài)上比初始模型稍大。

（3）Vyz反演結果物性值分布較集中，值也和初始模型最為接近，但基本上看不出初始模型的分布形態(tài)。

（4）Vzz的反演效果最差。

從下頁圖10分析得到，布格重力異常反演結果能夠大致表現(xiàn)出目標體的位置，但密度值略小且不集中。

圖11（見下頁）為全張量反演結果切片圖，可以看出全張量反演結果的密度值要略小，但仍較準確地體現(xiàn)出了異常體的密度，其幾何形態(tài)特征也與反演模型較接近，目標體所在位置也得到了較好的體現(xiàn)。因此可以認為，全張量反演綜合了各個張量分量反演的優(yōu)點。

4 結語

重力張量各個不同分量在反演過程中表現(xiàn)出對目標體不同敏感度，如目標體埋深、物性參數(shù)、目標形態(tài)等。而全張量反演能夠包含更多的場源信息，相對于單分量反演和布格重力異常反演，其反演結果有著更高的分辨率，且在識別場源特征方面有著更好的表現(xiàn)。

［1］L MARTIN.Advances and challenges in the development and deployment of gravity gradiometer systems［C］.EGM 2007International Workshop：Innovation in EM，Grav and Mag Methods：a new Perspective for Exploration.Capri，Italy，2007.

［2］MURPHY C A.Interpreting FTG gravity data using horizontal Tensor components［C］.EGM 2007International Workshop－Innovation in EM，Grav and Mag methods：new Perspective for Exploration，Ca－pri，Italy，2007.

［3］ALEXANDER DROUJININE，ALEXANDERR VASILEVSKY，RUSS EVANS.Feasibility of using full tensor gradient（FTG）data for detectionof local lateral density contrasts during reservoir monitoring［J］.Geophys，2007 （169）：795.

［4］XIONG LI，MICHEL CHOUTEAU.Three－Dimensional gravity modeling in all space［J］.Surveys in Geophysics，1998，19（4）：339.

［5］HINOJOSA J，K MICHUS.Hilbert transform of gravity gradient profiles：Special cases of the general gravity－gradient tensor in the Fourier transform domain［J］.Geophysics，2002：67.

［6］D W VASCO，C TAYLOR，Inversion of airborne gravity gradient data，southwestern Oklahoma［J］.Geophysics，1991，156（1）：90.

［7］MICHAEL S，ZHDANOV，ROBERT ELLIS，et al.Three－dimensional regularized focusing inversion of gravity gradient tensor component data［J］.Geophysics，2004，69（4）：925.

［8］JULIO CESAR SOARES DE OLIVEIRA LYRIO，LUIS TENORIO，YAOGUO LI，Efficient automatic denoising of gravity gradiometry data［J］.Geophysics，2004，69（3）：772.

［9］JAMES WHILE，ANDREW JACKSON，DIRK SMIT，et al.Spectral analysis of gravity gradiometry profiles［J］.Geophysics，2006，171（1）：11.

［10］KEVIN L.MICKUS，JUAN HOMERO HINOJOSA.The complete gravity gradient tensor derived from the vertical component of gravity：a Fourier transform technique［J］.Journal of Applied Geophysics：2001（46）：159.

［11］GARCIA－ABDESLEM J.Inversion of the power spectrum from gravity anomalies of prismatic bodies［J］.Geophysics，1995（60）：1698.

［12］CHANGYOU ZHANG， MARTIN F， MUSHAYANDEBVU，et al.Euler deconvolution of gravity tensor gradient data［J］.Geophysics，2000，65（2）：512.

［13］余欽范，樓海.水平梯度法提取重磁源邊界位置［J］.物探化探計算技術，1994，16（4）：363.

［14］王想，李桐林.Tilt梯度及其水平導數(shù)提取重磁源邊界位置［J］.地球物理學進展，2004，19（3）：625.

［15］魏偉，劉天佑.梯度法解釋復雜二維斷裂重力異常［J］.物探與化探，2005，29（8）：347.

［16］趙德軍，袁可佳.重力梯度張量計算重力場元［J］.2005，25（6）：12.

［17］李姍姍，吳曉平，吳星.重力矢量及張量信息在地球重力場中的應用［J］.測繪學院學報，2005，22（6）：94.

［18］劉尚余，周方.快速傅里葉變換在衛(wèi)星重力梯度測量中的應用［J］.華南理工大學學報：自然科學版，1997，25（7）：40.

［19］KENNDY J，EBERHART R C.Particle swarm optimization［C］.Proc IEEE Int Conf on Neural Networks，Perth，W A，Australia，1995.

［20］吳招才，劉天佑.板狀體磁異常數(shù)據(jù)反演的PSO算法［J］.物探與化探，2009，33（2）：194.

［21］邱寧，劉慶生，曾佐勛，等.基于混沌－粒子群優(yōu)化的磁法數(shù)據(jù)非線性反演方法［J］.地球物理學進展，2010，25（6）：2150.

［22］FEDI M，F(xiàn)ERRANTI L.Florio G，Giori I.Understanding the structural setting in the Southern Apennines（Italy）：insight from Gravity Gradient Tensor［J］.Tectonophysics，2005（397）：21.

［23］曾思紅.重力梯度張量正演研究及邊界提?。跠］.長沙：中南大學，2010.

［24］郭文斌.重力梯度張量的擬BP神經(jīng)網(wǎng)絡反演［D］.長沙：中南大學，2011.

［25］李剛毅，蔡涵鵬.基于粒子群優(yōu)化算法的波阻抗反演研究［J］.勘探地球物理進展，2008，31（3）：187.