基于Cokriging的三維重力張量隨機(jī)反演

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

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083)

?

基于Cokriging的三維重力張量隨機(jī)反演

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

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083)

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

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

0 引言

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

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

1 方法原理

1.1 重力張量模型正演

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

(1)

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

(2)

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

(3)

(4)

(5)

(6)

(7)

(8)

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

Vαβ=Sαβρ

(9)

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

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

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

(10)

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

CVVΛ=CVρ

(11)

最后利用一個(gè)理想的權(quán)系數(shù)可以從重力張量場(chǎng)中求解出密度的估計(jì)值：

ρ*=ΛTV

(12)

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

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

(13)

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

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

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

CVV=SCρρST+C0

(14)

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

CVρ=SCρρ

(15)

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

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

(16)

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

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

(17)

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

ρ*=ΛTV+ΦTρF

(18)

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

(19)

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

(20)

其中：z是網(wǎng)格單元的中心埋深；β一般情況下取一個(gè)接近2的數(shù)；z0為一個(gè)取決于長(zhǎng)方體單元尺寸和埋深的常數(shù)。通過調(diào)節(jié)z0和β的值，合理補(bǔ)償?shù)叵聵?gòu)造在深度上的分布。

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

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

(21)

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

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

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

(22)

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

(23)

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

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

(24)

2 模型算例

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

2.1 模型一

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

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

圖1 模型一圖

圖2 模型一重力張量正演結(jié)果

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

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

2.2 模型二

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

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

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

圖4 模型二圖

圖5 模型二重力張量正演結(jié)果

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

3 結(jié)論

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

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

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

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

[1] 王燦，朱自強(qiáng)，魯光銀，等. 重力梯度張量解析信號(hào)的正演模擬[J]. 中國(guó)有色金屬學(xué)報(bào),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ù)的三維反演方法與應(yīng)用[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] 陳少華. 基于預(yù)條件共軛梯度法的重力梯度張量反演研究[D]. 長(zhǎng)沙:中南大學(xué),2012. CHEN SH H. Inversion research of gravity gradient tensor based on preconditioned conjugate gradient.Changsha: Central South University,2012.(In Chinese)

[8] 宋慶磊，杜德文，丁明，等. 協(xié)同克里格方法在東海表面溫度場(chǎng)數(shù)據(jù)插值中的應(yīng)用[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)南，李維明. 實(shí)用地質(zhì)統(tǒng)計(jì)學(xué)[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]. 首都師范大學(xué)學(xué)報(bào):自然科學(xué)版,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]魯紅英，肖思和，何建軍，等. 三維地震資料儲(chǔ)層參數(shù)建模方法研究——分形協(xié)同克里格法[J]. 成都理工學(xué)院學(xué)報(bào),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)計(jì)學(xué)及其應(yīng)用[M].徐州:中國(guó)礦業(yè)大學(xué)出版社,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

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

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

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

1001-1749(2015)04-0428-09

P 631.1

A

10.3969/j.issn.1001-1749.2015.04.04