基于粒子群算法的重力張量反演研究

孫中科，朱自強(qiáng)，魯光銀，曹書錦

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

0 前言

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

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

1 正演基礎(chǔ)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 反演理論

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

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

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

3.2 目標(biāo)函數(shù)

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

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

3.3 模型試算

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

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

對該目標(biāo)體分別對各個張量單分量、布格重力異常和全張量進(jìn)行反演，反演耗時約460s～490s（因初始解是隨機(jī)設(shè)置的，同一張量反演的耗時也會因?yàn)槌跏冀庠O(shè)置的差別而略有不同），反演結(jié)果見圖2至圖5、及下頁圖6～圖9，并對比其效果。

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

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

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

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

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

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

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

4 結(jié)語

重力張量各個不同分量在反演過程中表現(xiàn)出對目標(biāo)體不同敏感度，如目標(biāo)體埋深、物性參數(shù)、目標(biāo)形態(tài)等。而全張量反演能夠包含更多的場源信息，相對于單分量反演和布格重力異常反演，其反演結(jié)果有著更高的分辨率，且在識別場源特征方面有著更好的表現(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］.物探化探計(jì)算技術(shù)，1994，16（4）：363.

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

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

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

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

［18］劉尚余，周方.快速傅里葉變換在衛(wèi)星重力梯度測量中的應(yīng)用［J］.華南理工大學(xué)學(xué)報：自然科學(xué)版，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］.地球物理學(xué)進(jìn)展，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］曾思紅.重力梯度張量正演研究及邊界提?。跠］.長沙：中南大學(xué)，2010.

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

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