重力梯度張量的預條件共軛梯度法反演

陳少華，朱自強，魯光銀，曹書錦

(中南大學 地球科學與信息物理學院，湖南 長沙，410083)

重力勘探方法在區(qū)域地質調查、成礦遠景區(qū)劃分評價、水文及工程地質勘察等方面發(fā)揮重要作用。常規(guī)重力勘探只觀測重力位的鉛垂一階導數(shù)，而重力梯度勘探則觀測重力位的二階導數(shù)，因此，相比傳統(tǒng)的重力測量，重力梯度測量有更高的分辨率[1]，更能夠清晰地反映場源體邊界形態(tài)，且其張量數(shù)據(jù)較為豐富，能夠提高異常體的解釋精度。作為一種前沿性的重力測量方法，國外學者不僅對重力梯度測量儀器進行了大量開發(fā)，還對采集的數(shù)據(jù)進行了大量研究。Zhang等[2]利用歐拉反褶積法解釋了重力全張量數(shù)據(jù)；Zhdanov等[3]利用重加權正則化共軛梯度法，采用聚焦正則化對目標函數(shù)進行約束，數(shù)值模擬表明該算法對尖銳邊界有較強的反演效果；Wilson等[4]針對巨大的航空重力梯度測量數(shù)據(jù)，采用并行算法，利用共軛梯度法對其進行反演，該算法在Vredefort得到了具體應用。此外，國外對重力張量數(shù)據(jù)的去噪、譜分析等方面的研究也取得了一定的成果[5-7]。在國內(nèi)，姚長利等[8]針對三維物性反演中存在的困難和問題，提出三維物性反演的隨機子域方法技術；郭良輝等[9]提出了基于重力梯度數(shù)據(jù)異常分離的三維相關成像方法提高成像分辨率,數(shù)值模擬結果表面該方法能夠使重力梯度數(shù)據(jù)具有良好的橫向及縱向分辨率。與國外相比，國內(nèi)對這一方面的研究較少，目前還沒有開發(fā)出重力梯度張量測量儀器，對重力梯度張量數(shù)據(jù)的解釋也較少。相較于非線性反演方法，線性反演算法出現(xiàn)的時間較長，應用較為廣泛。傳統(tǒng)線性反演算法中雅可比矩陣的存儲需占用大量內(nèi)存，計算該矩陣與其轉置相乘將耗費大量時間且占用大量內(nèi)存，而共軛梯度法只要計算兩次雅可比矩陣與向量的乘積[10-11]，計算時間將大量減少且所需內(nèi)存小，有利于大規(guī)模計算。針對雅可比矩陣的條件數(shù)較大而產(chǎn)生病態(tài)解的問題，首先對方程組進行預條件處理[12-13]，通過減小該矩陣的條件數(shù)來增強線性方程組的穩(wěn)定性。本文作者在目標函數(shù)中加入粗糙度函數(shù)及深度加權函數(shù)，利用預條件共軛梯度法對目標函數(shù)進行反演，數(shù)值模擬結果證明了反演的有效性；針對重力梯度張量采集的數(shù)據(jù)有5個分量是獨立的，充分利用重力梯度張量數(shù)據(jù)的豐富性，利用5個獨立分量進行反演，取得了較好的效果。

1 正演

首先，在笛卡兒坐標系下(圖1)，設地面上某點O為坐標原點。根據(jù)萬有引力定律，一個剩余密度為ρ、體積為V的地質體，其在外部空間任意一定產(chǎn)生的引力位為：

式中：r為觀測點 P(x,y,z)到異常點 Q(ξ,η,ζ)的矢量；G為萬有引力常量。

式中：g，gx，gy和gz均為笛卡兒坐標系(x,y,z)的函數(shù)；g為U(r)的一次導數(shù)；U(r)的二階導數(shù)可表示為：

由于在地球外部，位場滿足拉普拉斯方程Gxx+Gyy+Gzz=0，以及重力場的無旋性：Gxy=Gyx，Gxz=Gzx，Gyz=Gzy，因此，在重力梯度張量中只有5個分量是獨立的。

圖1 笛卡兒坐標系下的長方體模型及地下空間剖分Fig.1 Rectangular model and subsurface division in Cartesian coordinates

將探測區(qū)域的地下空間離散成一系列大小相等的長方體[14]，賦予每個小長方體1個固定的剩余密度，在探測區(qū)域建立笛卡兒坐標系(圖1)，x軸指向正北方向，y軸指向東方向，z軸垂直向下，則地面任意觀測點各梯度分量的異常值可表示為：

式中：iSαβ代表某重力梯度分量中第i個長方體的幾何架構，稱為位置函數(shù)；ρi為第i個長方體的剩余密度；M為所劃分長方體的個數(shù)。

則由地下空間離散的所有長方體在地面所有觀測點產(chǎn)生的重力梯度異常值滿足線性關系：

2 反演

2.1 粗糙度約束

重力梯度張量的反演便是求解線性方程組(5)。由于采集的數(shù)據(jù)向量d遠少于要反演的未知數(shù)m，因此該反演是個欠定問題，方程組的解不唯一且不穩(wěn)定，必須對目標函數(shù)加上適當?shù)募s束來縮小解的范圍。本文中通過引入粗糙度函數(shù)[15]，來對模型空間進行約束，目標函數(shù)可表示為：

定義模型粗糙度R為m(r)在x，y和z方向一階偏微分的平方和，即

2.2 深度加權

式中：z為塊體單元中心點埋深；z0和β為常數(shù)。β一般情況下取 2～4；z0取決于塊體單元的尺寸以及觀測面的高度。

將式(8)及最小化模型構造代入式(7)：

式中：αs，αx，αy和 αz為模型目標函數(shù)中各項的加權因子。

對式(9)按照模型區(qū)域的網(wǎng)格剖分進行離散化，用有限差分代替微分算子，寫成如下矩陣形式：

每一個矩陣分量都可以寫成2個獨立的矩陣與1個系數(shù)矩陣的乘積，即

式中：Ri為有限差分算子；Z為離散化的深度加權對角線矩陣。將式(10)代入式(6)：

式(12)可寫成方程組：

式(13)可寫成：

2.3 預條件共軛梯度法

由于模型參數(shù)的規(guī)模較大，在實際計算中，收斂的次數(shù)越少越好，而方程的收斂速度由雅可比矩陣的條件數(shù)決定的。位場中雅可比矩陣的條件數(shù)通常巨大，同時，van Decar等[18]指出在目標函數(shù)中增加粗糙度約束將減緩收斂速度，為了避免過多的收斂次數(shù)及增強解的穩(wěn)定性，將常規(guī)的方程AtAm=Atb改進成：

因此，預條件共軛梯度法的計算步驟：

(1) 給定初始值ε，令m0=0；k=0；r0=Atb；

(4) 計算 rk=rk-1-αkAtqk，若||rk||2＜ε，則停止計算，否則返回第(2)步。

其中：S為預條件因子。若 S≈(AtA)-1，則 SAtA≈I(其中，I為單位矩陣)，那么預條件因子S就可以大大改善式(14)的條件數(shù)，使得特征值集中分布在對角線附近，有效地提高迭代速度，這便是預條件共軛梯度法。但是，這樣求解S是非常困難的，在實際應用中，用簡單的深度加權函數(shù)即式(8)作為預條件因子[12]，能很好地改善式(14)的條件數(shù)，實際模型反演結果證明了這一點。

3 數(shù)值模擬

3.1 單分量反演

將場源空間劃分為20×20×10=4 000個單元格，每個單元格沿x，y和z軸方向的長度均為50 m。如圖2所示，地下空間存在2個異常體，2個異常體的長×寬×高均為200 m×200 m×200 m，異常體頂板埋深均為125 m，底板為325 m，異常體的剩余密度為1.0 t/m3?？偣灿?0×20=400個采集數(shù)據(jù)，數(shù)據(jù)采集點位于表層，沿x和y軸間距都為50 m。

圖 2 正演模型和 Gxx，Gxy，Gxz，Gyy，Gyz和Gzz的反演結果Fig.2 Forward model and inversion results of Gxx, Gxy, Gxz, Gyy, Gyz and Gzz

定義絕對擬合差：

反演參數(shù)的設置中，取 αx=αy=αz=1，αx=0.00 005，z0=20，β=4，收斂終止的條件是上一次迭代的絕對擬合差大于后一次的絕對擬合差。從以上各個單分量反演結果可知：每個分量基本能夠反演出異常地質體的位置和形態(tài)，反演出來的剩余密度也都不同程度地接近真實值，但都出現(xiàn)了大小不等的負值，其中以 Gxx和 Gxz出現(xiàn)的負值較大，而 Gxy，Gyy，Gzz和 Gyz反演出來的負值相對較??；同時可知Gxx，Gxz，Gxy和 Gyy分量的反演結果較為發(fā)散，在深度方向上的延伸較大，不能完全反演出異常點的位置；而Gyz及Gzz分量的反演結果集中于異常體周圍，且準確地反映了目標體的深度。

為了各個分量的比較，定義相對擬合差Rms：

圖3所示為各個分量在不同迭代次數(shù)下的相對擬合差。從圖3可知：Gxz分量的反演非常不穩(wěn)定，隨著迭代的進行Rms上下波動比較厲害，Gxx分量也出現(xiàn)了小范圍的上下波動，這可能是以上2個分量反演結果中出現(xiàn)較大負值的原因；Gxy及Gyy分量基本上隨著迭代的進行 Rms逐漸變小，但減小的速度較慢，最終其Rms仍然較大；Gyz及Gzz分量不僅隨著迭代的進行Rms快速減小，而且最終的Rms也比較小。

圖3 各個分量在不同迭代次數(shù)下的相對擬合差Fig.3 Relative misfit of all single component in different numbers of iteration

3.2 全張量反演

雖然理論上可以證明同一方法不同分量并不能改變由位場性質決定的多解性[19]，這由于引起位場多解性的一個重要原因是數(shù)據(jù)的不完整和不精確，而觀測得到的數(shù)據(jù)總是受限于一定網(wǎng)格距和一定的精度。盡管如此，同一高度不同分量的觀測數(shù)據(jù)在完整性上仍有所補充，并減少了誤差引起的多解性。重力梯度張量相比于重力勘探，不僅有更高的分辨率，還在于其更為豐富的數(shù)據(jù)量，充分利用5個獨立分量的信息將有利于提高反演的精度。

本文采用5個獨立分量Gxx，Gxy，Gxz，Gyy及Gzz進行聯(lián)合反演，并對反演數(shù)據(jù)加入10%的噪聲干擾，檢驗反演的抗噪性能。取 αx=αy=αz=1，αs=0.000 05，z0=10，β=4，雖然隨著迭代的進行，數(shù)據(jù)的絕對擬合差能不斷地降低，但可能導致過擬合而引起多余構造，故本次收斂中取最多的迭代次數(shù)為30次，觀測數(shù)據(jù)的正演模型見圖2(a)。

與沒有加入高斯噪聲的圖 2中各個單分量相比較：圖4(a)所示的反演的異常體不僅集中，而且和異常體的真實深度吻合，驗證了同一種方法的不同分量的聯(lián)合反演可以減少多解性，提高解的精度。圖4(b)中加入 10%的高斯噪聲后反演的異常形態(tài)仍較為明顯，只是在異常體的周圍出現(xiàn)一些小異常，但可以忽略不計，表面5個分量聯(lián)合反演有較好的抗噪性。

3.3 合成模型反演

現(xiàn)設計1個Y型巖脈，該模型曾被姚長利等[8-9,20]應用于位場相關領域的分析。圖5(a)和(b)所示分別為深度120 m的深度切片和圖5(a)中A-A′的斷面切片，圖5(c)所示為120 m深度剖面的反演結果，圖5(d)所示為圖5(a)中A-A′斷面切片的反演結果，白色方框為真實異常體的位置。異常區(qū)Ⅰ(顏色較淺)的異常體其剩余密度為0.8 t/m3，深度方向延伸不大；異常區(qū)Ⅱ(顏色較深)的異常體其剩余密度為1.0 t/m3，深度方向延伸大，兩巖脈的延伸方向相反。將地下空間剖分為20×20×10=4 000個長方體，每個長方體的大小都為50 m×50 m×50 m。正演計算的觀測點位于地表處，共有20×20=400個觀測點，在觀測數(shù)據(jù)中加入3%的高斯噪聲，采用5個獨立分量對地下異常體進行聯(lián)合反演。

反演參數(shù)的設置中，取 αx=αy=αz=1，αx=0.000 05，z0=20，β=4，最多迭代次數(shù)為30次，加入先驗參數(shù)范圍約束，當解小于0時令其為0，當解大于1.1時令其為1.1。整個反演過程的迭代次數(shù)為30次，消耗的時間為21.742 587 s，相對擬合差為0.026 371。

圖 5(c)所示較好地反映了異常體的形態(tài)和位置，左邊巖脈的形態(tài)和位置基本吻合，右邊巖脈沿y軸方向完全吻合，在x軸方向稍微擴大了一點，但基本上與異常體的真實形態(tài)一致；圖5(d)中基本重構出了模型的異常形態(tài)，其異常走向與Y型巖脈一致；結合圖5(c)和5(d)可發(fā)現(xiàn)：反演算法在淺部取得了較清晰的結果，而在深部反演的分辨率則有所下降，表明雖然基于深度加權的重力梯度張量反演能夠使物性的參數(shù)值大致能夠分布在相應的位置，但重構模型的分辨率不可避免地隨著深度的增加而降低。以上結果表明：基于粗糙度約束和深度加權的重力梯度張量反演，能夠在一定程度上較好地重構模型的異常形態(tài)。

4 結論

(1) 提出了基于預條件共軛梯度法的重力梯度張量反演，收斂速度快，耗時較短，可用于數(shù)據(jù)量較大的反演計算。

(2) 相對于單分量反演，充分利用 5個獨立分量的信息進行聯(lián)合反演，能夠較準確地反演出異常體的形態(tài)及剩余密度特征。

(3) 僅用理想模型數(shù)據(jù)對本文的算法進行了驗證，需要進一步采用現(xiàn)場采集的數(shù)據(jù)進行驗證。

[1] 曾華霖. 重力場與重力勘探[M]. 北京: 地質出版社, 2005:68-71.ZENG Hualin. gravitational field and gravitational prospecting[M]. Beijing: Geological Publishing House, 2005:68-71.

[2] Zhang C, Mushayandebvu M F, Reid A, et al. Euler deconvolution of gravity tensor gradient data[J]. Geophysics,2000, 65(2): 512-520.

[3] Zhdanov M S, Ellis R, Mukherjee S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2004, 69(4): 925-937.

[4] Wilson G, Cuma M, Zhdanov M S. Massively parallel 3D inversion of gravity and gravity gradiometry data[J]. Preview,2011, 152(6): 29-34.

[5] Lyrio L S, Tenorio L, Li Y G. Efficient automatic denoising of gravity gradiometry data[J]. Geophysics, 2004, 69(3): 772-782.

[6] While J, Jackson A, Smit D, et al. Spectral analysis of gravity gradiometry profiles[J]. Geophysics, 2006, 71(1): J11-J22.

[7] Afif H. Understanding gravity gradients: Atutorial[J]. The Leading Edge, 2006, 25(8): 942-947.

[8] 姚長利, 鄭滿元, 張聿文. 重磁異常三維物性反演隨機子域法方法技術[J]. 地球物理學報, 2007, 50(5): 1576-1583.YAO Changli, ZHENG Manyuan, ZHANG Yuwen. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces[J]. Chinese Journal of Geophysics, 2007, 50(5):1576-1583.

[9] 郭良輝, 孟曉紅, 石磊, 等. 重力和重力梯度數(shù)據(jù)三維相關成像[J]. 地球物理學報, 2009, 52(4): 1098-1106.GUO Lianghui, MENG Xiaohong, SHI Lei, et al. 3-D correlation imaging for gravity and gravity gradiometry data[J].Chinese Journal of Geophysics, 2009, 52(4): 1098-1106.

[10] 吳小平, 汪彤彤. 利用共軛梯度方法的電阻率三維反演的若干問題研究[J]. 地震地質, 2001, 23(2): 321-327.WU Xiaoping, WANG Tongtong. Study on some problems for 3D resistivity inversion using conjugate gradient[J]. Seismology and Geology, 2001, 23(2): 321-327.

[11] 吳小平, 徐國明. 利用共軛梯度法的電阻率三維反演研究[J].地球物理學報, 2000, 43(3): 420-427.WU Xiaoping, XU Guoming. Study on 3-D resistivity inversion using conjugate gradient method[J]. Chinese Journal of Geophysics, 2000, 43(3): 420-427.

[12] Pilkington M. 3-D magnetic imaging using conjugate gradient[J].Geophysics, 1997, 62(4): 1132-1142.

[13] Pilkington M. 3D magnetic data-space inversion with sparseness constraints[J]. Geophysics, 2009, 74(1): L7-L15.

[14] 駱遙, 姚長利, 薛典軍, 等. 2.5D地質體重磁異常無解析奇點正演計算研究[J]. 石油地球物理勘探, 2009, 44(4): 487-493.LUO Yao, YAO Changli, XUE Dianjun, et al. Analytic singularity-free forward computation of gravity and magnetic anomaly for 2.5-D geologic body[J]. Oil Geophsical Prospecting,2009, 44(4): 487-493.

[15] Constable S C, Parker R L, Constable C G. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[J]. Geophysics, 1987, 52(3):289-300.

[16] Li Y G, Oldenburg D W. 3-D inversion of magnetic data[J].Geophysics, 1996, 61(2): 394-408.

[17] 姚姚. 地球物理反演基本理論與應用方法[M]. 武漢: 中國地質大學出版社, 2002: 57-59.YAO Yao. Basic theory and application of geophysical inversion methods[M]. Wuhan: China University of Geosciences Press,2002: 57-59.

[18] van Decar J, Snieder R. Obtaining smooth solution to large,linear, inverse problems[J]. Geophysics, 1994, 59(5): 818-829.

[19] 管志寧. 地磁場與磁力勘探[M]. 北京: 地質出版社, 2005:199-204.GUAN Zhining. Magnetic field and magnetic exploration[M].Beijing: Geological Publishing House, 2005: 199-204.

[20] Li Y G, Oldenburg D W. 3-D inversion of gravity data[J].Geophysics, 1998, 63: 103-119.