黃天統(tǒng)， 朱自強， 魯光銀， 曹書錦

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

?

基于Cokriging的三維重力張量隨機反演

黃天統(tǒng)， 朱自強*， 魯光銀， 曹書錦

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

重力張量是重力位的二階導數(shù)，與傳統(tǒng)的重力測量相比,重力張量具有更高的分辨率。協(xié)同克里格法是地質(zhì)統(tǒng)計學中常見的一種方法，使用協(xié)同克里格法對三維重力張量各分量數(shù)據(jù)進行隨機反演。由于位場數(shù)據(jù)在深度方向上的分辨率較低，因此需通過使用深度加權函數(shù)和鉆孔數(shù)據(jù)來解決這一問題，對比了帶深度加權和帶鉆孔數(shù)據(jù)的協(xié)同克里格反演結果。反演結果表明，協(xié)同克里格法能夠較好地反演重力張量各分量數(shù)據(jù)，兩種方式都能夠較好地反演出地下異常位置，帶深度加權的協(xié)同克里格法反演出來的密度值與實際值有一定的偏差，而帶鉆孔數(shù)據(jù)的協(xié)同克里格法能夠較好地反演出剩余密度值，并且對于不同埋深的異常體都能很好地反映。

協(xié)同克里格； 重力張量； 反演； 鉆孔數(shù)據(jù)； 深度加權函數(shù)

0 引言

作為一種前沿性的重力測量方法，重力張量測量在國外開展已有十余年。與常規(guī)的重力測量相比，重力張量測量的是重力位的二階導數(shù)，對地下介質(zhì)密度的變化反映更加靈敏，能夠更加直接地突出目標體的邊界[1]，因此選擇一種合理的數(shù)值算法對重力梯度張量進行反演具有重要意義。LI Xiong[2]提出的全空間的重力梯度張量計算公式，為重力梯度張量理論計算提供了依據(jù)，也是重力梯度張量反演理論的基礎；李耀國[3-4]在重磁數(shù)據(jù)三維反演時，為了提高反演在深度方向的分辨率，引入了深度加權函數(shù)；Zhdanov等人[5]實現(xiàn)了三維重力梯度張量數(shù)據(jù)聚焦正則化反演，反演結果對異常體的尖銳邊界反映明顯；王浩然[6]將李耀國的方法應用到重力梯度張量反演，構建了聯(lián)合反演目標函，并對反演密度進行約束，反演結果能夠很好地反映地下異常體位置和剩余密度；陳少華[7]利用預條件共軛梯度對重力梯度張量各分量數(shù)據(jù)進行聯(lián)合反演取得了很好的效果。

地質(zhì)統(tǒng)計學是數(shù)學地質(zhì)領域中發(fā)展迅速且應用廣泛的一門學科，它是在區(qū)域化變量理論的基礎上，運用變差函數(shù)這一主要工具，研究變量在空間分布的隨機性與結構性[8-9]。協(xié)同克里格法[10]則把區(qū)域化變量的最佳估值方法從單一屬性發(fā)展到兩個以上的協(xié)同區(qū)域化屬性。國外有部分學者將地質(zhì)統(tǒng)計學方法運用到地球物理反演方面，取得了一定的成果。Bosch[11-12]對密度的后驗概率方程應用蒙特卡羅方法較好地反演出初始模型；Chasseriau[13]應用先驗模型的協(xié)方差對重力數(shù)據(jù)進行了反演。協(xié)同克里格法能夠簡便地結合鉆孔數(shù)據(jù)，使得反演結果更加準確可靠。Shamsipour[14-15]等人利用地質(zhì)統(tǒng)計學結合鉆孔數(shù)據(jù)對重力數(shù)據(jù)進行了三維cokring反演和條件模擬，反演結果能夠準確地確定異常邊界和剩余密度。Erwan Gloaguen等人[16]將協(xié)同克里格法運用于鉆孔地質(zhì)雷達波速反演取得了很好的效果。作者將協(xié)同克里格方法應用于重力張量的反演，并結合深度加權和鉆孔數(shù)據(jù)，反演結果能夠較好地反演出地下異常體的邊界和剩余密度值。

1 方法原理

1.1 重力張量模型正演

重力張量是重力位的二階導數(shù)，在笛卡爾坐標系下，根據(jù)牛頓萬有引力定律，一個剩余密度為ρ，體積為V的地質(zhì)體，在它外部空間中的任意一點的重力位為：

(1)

其中：G是萬有引力常量；r是任意觀測點P(x,y,z)到異常點Q(ξ,η,ζ)的方向矢量。那么重力張量表達式寫成矩陣形式：

(2)

式中：Γ表示全重力梯度張量。重力張量共有9個分量，由于在無源空間里重力的旋度和散度都為"0"，因此重力梯度張量中只有5個獨立分量Vxx、Vxy、Vxz、Vyz、Vzz。對于長方體，重力張量各分量的離散正演公式[2]為：

(3)

(4)

(5)

(6)

(7)

(8)

若地下剖分的長方體單元格的數(shù)量為m，測點數(shù)為n，將重力張量各分量寫成統(tǒng)一的表達式：

Vαβ=Sαβρ

(9)

其中：α、β分別取值為x、y、z，不同的取值對應不同的重力張量分量；Vαβ為n×1維的觀測數(shù)據(jù)向量；ρ為m×1維單元格剩余密度的值；Sαβ為n×m維的雅克比矩陣；矩陣S的元素Sij表示第j個單元格在第i個測點的響應。

1.2 協(xié)同克里格原理

協(xié)同克里格方法是地質(zhì)統(tǒng)計學中重要而且十分有效的方法之一，它既能反映統(tǒng)計相關信息，又能反映空間位置的相關信息[17]。協(xié)同克里格方法是通過運用主要變量和次要變量的空間相關性來提高主變量的估計。在這里剩余密度ρ以及其估計值ρ*是主要變量，重力張量V是次要變量。假設重力張量和剩余密度的期望是為“0”，即E[ρ]=E(V)=0，若期望不為零，可先減去平均值，最后再補償?？梢詮?式(10)矩陣的對角元素上求出估計方差：

(10)

其中：CVV是重力張量的協(xié)方差矩陣；Cρρ是密度的協(xié)方差矩陣；CρV是密度和重力張量的交互協(xié)方差矩陣；Λ是權系數(shù)矩陣；ρ和V是多維的隨機變量。最小化式(10)的估計方差可以得到一個簡單的協(xié)同克里格的解：

CVVΛ=CVρ

(11)

最后利用一個理想的權系數(shù)可以從重力張量場中求解出密度的估計值：

ρ*=ΛTV

(12)

協(xié)同克里格方差就可以通過公式(13)求解：

σck=diag(Cρρ-ΛTCVρ)

(13)

矩陣A=Cρρ-ΛTCVρ中的非對角元素表示估計誤差的協(xié)方差。

1.3 協(xié)同克里格反演

從公式(9)可以看出，剩余密度和重力張量的協(xié)方差矩陣是線性相關的：

CVV=SCρρST+C0

(14)

其中:C0是重力張量觀測誤差的協(xié)方差矩陣，一般情況下，C0是一個對角陣。

CVρ=SCρρ

(15)

如果不考慮觀測誤差，即C0=0，就有式(16)成立:

=S(SCρρ)T(SCρρST)-1V=V

(16)

從式(16)可以看出，在沒有誤差的情況下，利用協(xié)同克里格計算出來的剩余密度就等于觀測場的剩余密度。

由于重力張量數(shù)據(jù)和其他位場數(shù)據(jù)一樣，在深度方向上的分辨率較低，反演出來的構造都會集中在地表。因此要想加強深度方向上的分辨率還得進行一定的處理。如果已知鉆孔的數(shù)據(jù)資料，則可以提高深度方向上的分辨率，如果鉆孔處的剩余密度值為ρF，通過擴展式(11)可以得到類似的公式：

(17)

CVρF為重力張量與鉆孔剩余密度的交互協(xié)方差，矩陣CρFV=CVρFT為鉆孔剩余密度與重力張量的交互協(xié)方差，CρFρ為鉆孔剩余密度與區(qū)域剩余密度的交互協(xié)方差，Λ和Φ是最小化方差估計的權系數(shù)，Giroux[18]推導了式(17)中各變量的求取方法。

ρ*=ΛTV+ΦTρF

(18)

那么協(xié)同克里格方差就可以通過式(19)求解：

(19)

在沒有鉆孔數(shù)據(jù)的情況下，則可以通過引入深度加權因子[4]來補償深度方向上的衰減,從而增強深度方向的分辨率。

(20)

其中：z是網(wǎng)格單元的中心埋深；β一般情況下取一個接近2的數(shù)；z0為一個取決于長方體單元尺寸和埋深的常數(shù)。通過調(diào)節(jié)z0和β的值，合理補償?shù)叵聵嬙煸谏疃壬系姆植肌?/p>

1.4 模型協(xié)方差估計

假設E[ρ]=E(V)=0，那么觀測的數(shù)據(jù)的協(xié)方差矩陣可以寫成[14]：

(21)

由于測區(qū)密度在水平方向被截斷，使協(xié)方差矩陣CVV、CVρ都是不穩(wěn)定的，因此傳統(tǒng)的估計方法不能直接運用。Asli等人[19]運用理論變差和實際變差圖像的方法(V-V圖像法)解決了這一問題。

1)假定一個剩余密度協(xié)方差矩陣Cρρ的初始模型。變差函數(shù)與協(xié)方差矩陣之間的關系式為：

C(h)=C(0)-γ(h)

(22)

γ(h)為模型的變差模型，例如式(23)為球體變差模型，h為空間距離，通過改變變差模型可以調(diào)節(jié)Cρρ。對于簡單的模型球體變差模型能夠很好地擬合，若當球體變差模型無法滿足擬合精度要求時，可以選擇其他變差模型[20]。

(23)

其中：C1為塊金常數(shù)；C1+C為基臺值；C為拱高；a為變程。

2)根據(jù)式(21)計算CVV的實驗值CVV-exp，根據(jù)式(14)計算CVV的理論值CVV-th，并將兩個矩陣重排成兩個向量νexp和νth。

(24)

2 模型算例

分別以單個長方體和水平方向上不同埋深、不同密度的兩個長方體為研究模型，利用協(xié)同克里格的方法對兩個模型進行反演分析。

2.1 模型一

將場源空間劃分為20×20×20=8 000個單元格，在x、y、z三個方向單元格的長度均為50 m，如圖1(a)所示，場源空間存在一個400 m×400 m×200 m的長方體，長方體頂部埋深為200 m，底部埋深為400 m，長方體水平中心位置為(475 m,475 m)，剩余密度為1 000 kg/m3,x、y方向點距為50 m，測線長度為1 000 m，鉆孔位置在(475 m,475 m)，鉆孔數(shù)據(jù)如圖1(b)所示，模型一重力張量正演結果如圖2所示。

運用V-V圖像法，取步長Nlag=100，調(diào)節(jié)模型協(xié)方差矩陣使理論協(xié)方差矩陣和實驗協(xié)方差相近，以Vzz的反演為例，通過擬合變差模型，最小化式(24)，得到球體變差模型C=7 000 (kg/m3)2，C1=0，ax=350 m，ay=350 m，az=350 m。反演結果如圖3所示。

圖1 模型一圖

圖2 模型一重力張量正演結果

圖3 模型一協(xié)同克里格反演結果(切片在y=400 m處)

對比圖3(a)和圖3(b)，對于單個的長方體，帶深度加權的協(xié)同克里格反演和帶鉆孔數(shù)據(jù)的協(xié)同克里格反演，都能夠比較準確地反演出地下異常的位置，但帶鉆孔數(shù)據(jù)的協(xié)同克里格反演對異常位置確定的效果，明顯優(yōu)于帶深度加權的協(xié)同克里格反演的結果。對于反演的密度值，帶深度加權的協(xié)同克里格反演出來的密度值與初始模型有較大的偏差，而帶鉆孔數(shù)據(jù)的協(xié)同克里格反演能夠較好地反演出密度值。觀察圖3(b)至圖3(f)可以看出，結合鉆孔數(shù)據(jù)，重力張量各分量反演出來的結果，都能較好地確定地下異常體的位置和剩余密度值。

2.2 模型二

如圖4(a)所示，地下空間水平方向存在兩個長方體，大小都為200 m×200 m×200 m，左側長方體水平中心位置為(275 m,500 m)，頂部埋深為200 m，底部埋深400 m，剩余密度為1 000 kg/m3；右側長方體水平中心位置為(675 m,500 m)，頂部埋深為100 m，底部埋深300 m，剩余密度為500 kg/m3，鉆孔位置在(300 m,500 m)和(700 m,500 m)，鉆孔數(shù)據(jù)如圖4(b)所示，模型二重力張量正演結果如圖5所示。

以Vzz的反演為例，通過擬合變差模型，最小化式(24)，得到球體變差模型C=5 000 (kg/m3)2，C1=0，ax=300 m，ay=300 m，az=400 m。反演結果如圖6所示。

對比圖6(a)和圖6(b)，當反演在水平方向不同埋深、不同密度的兩個長方體時，兩種方法都能夠很好地反演出地下異常的位置，對于反演的密度值，帶深度加權的協(xié)方差方法反演出的埋深較淺的異常體密度值比埋深較深的異常體密度值大，與初始模型相反，這是由于在模型二的情況下，埋深較深的異常體在地表產(chǎn)生的重力張量異常要比淺部異常體在地表產(chǎn)生重力張量值小。帶鉆孔數(shù)據(jù)的協(xié)同克里格反演能夠準確地確定兩個長方體的具體位置和密度值。觀察圖6(b)至圖6(f)可以看出，結合鉆孔數(shù)據(jù)，重力張量各分量反演出來的結果都能較好地確定地下異常體的位置及剩余密度值。

圖4 模型二圖

圖5 模型二重力張量正演結果

圖6 模型二協(xié)同克里格反演結果(切片在y=450 m處)

3 結論

作者提出了基于協(xié)同克里格法重力張量隨機反演，并通過長方體模型驗算了協(xié)同克里格的反演效果，并得出以下結論：

1)利用協(xié)同克里格法能夠反演重力張量各分量數(shù)據(jù)，反演結果能夠較好地確定地下異常體的位置。

2)利用帶深度加權函數(shù)的協(xié)同克里格方法反演出來的密度值跟初始模型有偏差，對多個異常體反演時，由于淺部異常的響應，帶深度加權函數(shù)的協(xié)同克里格方法反演出的較深部信息不明顯。

3) 利用鉆孔數(shù)據(jù)反演出來的結果能夠較好的確定剩余密度值，并且對不同埋深的異常體都能很好地反映。

[1] 王燦，朱自強，魯光銀，等. 重力梯度張量解析信號的正演模擬[J]. 中國有色金屬學報,2013,23(9):2479-2497. WANG C, ZHU Z Q,LU G Y,et al. Forward simulation of analytic signals of gravity gradient tensor. The Chinese Journal of Nonferrous Metals, 2013,23(9): 2479-2497.(In Chinese)

[2] LI X, CHOUTEAU M. Three-dimensional gravity modeling in all space[J]. Surveys in Geophysics,1998,19: 339-368.

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

[4] LI Y G, OLDENBURG D W. 3-D inversion of gravity data[J].Geophysics,1998,63(1):109-119.

[5] ZHDANOV M, ELLIS R, MUKHERJEE S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics,2004, 69(4):925-937.

[6] 王浩然，陳超，杜勁松.重力梯度張量數(shù)據(jù)的三維反演方法與應用[J]. 石油地球物理勘探,2013,48(3):474-481. WANG H R,CHEN CH,DU J S. 3-D inversion of gravity gradient tensor data and its applications.OGP,2013,48(3): 474-481.(In Chinese)

[7] 陳少華. 基于預條件共軛梯度法的重力梯度張量反演研究[D]. 長沙:中南大學,2012. CHEN SH H. Inversion research of gravity gradient tensor based on preconditioned conjugate gradient.Changsha: Central South University,2012.(In Chinese)

[8] 宋慶磊，杜德文，丁明，等. 協(xié)同克里格方法在東海表面溫度場數(shù)據(jù)插值中的應用[J]. 海岸工程,2011,30(3): 49-55. SONG Q L,DU D W,DING M,et al. Application of co-kriging method to interpolation of surface temperature data in East China sea .Coastal Engineering, 2011,30(3):49-55.(In Chinese)

[9] 侯景儒，尹鎮(zhèn)南，李維明. 實用地質(zhì)統(tǒng)計學[M]．北京：地質(zhì)出版社，1998. HOU J R,YIN ZH N,LI W M. Practical geostatistics.Beijing: Geological Publishing House,1998.(In Chinese)

[10]季青，余明. 基于協(xié)同克里格插值法的年均溫空間插值的參數(shù)選擇研究[J]. 首都師范大學學報:自然科學版,2010,31(4):81-87. JI Q,YU M. Study on parameters setting of ordinary cokriging interpretation to average annual temperature.Journal of Capital Normal University: Natural Science Edition,2010,31(4): 81-87.(In Chinese)

[11]BOSCH M, MCGAUGHEY J. Joint inversion of gravity and magnetic data under lithologic constraints[J].The leading edge,2001,20:877-881.

[12]BOSCH M, MEZA R, JIMéNEZ R, et al. Joint gravity and magnetic inversion in 3D using Monte Carlo methods[J]. Geophysics, 2006,71(4):153-156.

[13]CHASSERIAU P, CHOUTEAU M. 3D gravity inversion using a model of parameter covariance[J]. Journal of Applied Geophysics, 2003, 52:59-74.

[14]SHAMSIPOUR P, CHOUTEAU M,MARCOTTE D. 3D stochastic inversion of magnetic data[J]. Journal of Applied Geophysics,2011(73): 336-347.

[15]SHAMSIPOUR P, MARCOTTE D, CHOUTEAU M. 3D stochastic inversion of gravity data using cokriging and cosimulation[J]. Geophysics,2010,75(1):I1-I10.

[16]GLOAGUEN E, MARCOTTE E, GIROUX B ET AL. Stochastic borehole radar velocity and attenuation tomographies using cokriging and cosimulation[J]. Journal of Applied Geophysics,2007,62:141- 157.

[17]魯紅英，肖思和，何建軍，等. 三維地震資料儲層參數(shù)建模方法研究——分形協(xié)同克里格法[J]. 成都理工學院學報,2002,29(6):694-697. LU H Y,XIAO S H,HE J J,et al. Research of method on reservoir parameters model using 3D seismic data containing rock-fractal cokriging.Journal of Cheng Du University of Technology,2002,29(6):694-697.(In Chinese)

[18]GIROUX B, GLOAGUEN E, CHOUTEAU M. bh_tomo—a Matlab borehole georadar 2D tomography package[J]. Computers and Geosciences,2007,33:126-137.

[19]ASLI M, MARCOTTE D,CHOUTEAU M. Direct inversion of gravity data by cokriging[C].In: Kleingeld, W., Krige, D. (Eds.), 6th International Geostatistics Congress, Cape town. South Africa.2000:64-73.

[20]孫洪泉.地質(zhì)統(tǒng)計學及其應用[M].徐州:中國礦業(yè)大學出版社,1990. SUN H Q. Geostatistics and its application .Xuzhou:China University of Mining and Technology press,1990.(In Chinese)

3D stochastic inversion of gravity tensor based on Cokriging

HUANG Tian-tong, ZHU Zi-qiang*, LU Guang-yin, CAO Shu-jin

(School of Geosciences and Info-Physics, Central South University，Changsha 410083,China)

Gravity tensor is the second derivative of gravitational field. Compared with traditional gravitational exploration, gravity gradient tensor has a better resolution. Cokriging is a common method in geostatistics. In this paper, we applied cokriging to 3D stochastic inversion of each gravity tensor componentdata. Because potential field data have little resolution along the depth direction, we use depth weighting function and borehole data to solve this problem, and we compare the inversion results between both above methods. The inversion results show that the cokriging method can invert each gravity tensor component data. And both methods are capable of inverting the position of causative in the subsurface. In addition inversion results with depth weighting function cannot precisely recovered the anomalous bodies. But incorporating borehole date into inversion can improve the results with better resolution along the depth.

cokriging; gravity tensor; inversion; borehole date; depth weighting function

2014-07-01 改回日期：2014-08-14

國家自然科學基金(41174061，41374120); 中南大學自由探索計劃(2011QNZT011)

黃天統(tǒng)(1991-)，男，碩士，從事重磁相關領域的研究，E-mail:huang3200101@163.com。

*通信作者：朱自強(1964-)，教授，博士，主要從事重磁相關領域的研究,E-mail：13507319431@139.com。

1001-1749(2015)04-0428-09

P 631.1

A

10.3969/j.issn.1001-1749.2015.04.04