耿美霞， 黃大年， 于平*， 楊慶節(jié)

1 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢　430074 2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春　130026 3 中國(guó)科學(xué)院測(cè)量與地球物理研究所， 武漢　430077

基于協(xié)同克里金方法的重力梯度全張量三維約束反演

耿美霞1,2， 黃大年2， 于平2*， 楊慶節(jié)3

1 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢430074 2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春130026 3 中國(guó)科學(xué)院測(cè)量與地球物理研究所， 武漢430077

摘要隨著重力梯度全張量測(cè)量技術(shù)的日趨成熟，重力梯度全張量數(shù)據(jù)的三維反演技術(shù)日益受到重視與關(guān)注.全張量數(shù)據(jù)反演與重力數(shù)據(jù)反演一樣仍然面臨著嚴(yán)重的多解性問(wèn)題.本文將基于地質(zhì)統(tǒng)計(jì)學(xué)的協(xié)同克里金方法應(yīng)用于重力梯度全張量數(shù)據(jù)三維反演，建立了密度約束下的多變量協(xié)同克里金聯(lián)合反演方程，以降低反演的多解性.模型試驗(yàn)表明密度信息的加入能夠有效降低反演的多解性，提高反演結(jié)果的分辨率，尤其是縱向分辨率能夠得到顯著提高.最后對(duì)美國(guó)德克薩斯州一個(gè)巖丘區(qū)所獲得實(shí)際資料的應(yīng)用表明了本文方法的實(shí)用性.

關(guān)鍵詞重力梯度全張量； 三維約束反演； 協(xié)同克里金方法； 地質(zhì)統(tǒng)計(jì)學(xué)

1引言

重力梯度測(cè)量能夠觀測(cè)重力位的二階導(dǎo)數(shù)，因而，與傳統(tǒng)的重力測(cè)量相比，具有更高的分辨率，能夠更好地反映淺部異常的邊界(曾華霖，1999).在移動(dòng)平臺(tái)中(例如船和飛機(jī))，梯度儀不像重力儀那樣易受共模噪聲和大的加速的影響，因此與重力數(shù)據(jù)相比，重力梯度數(shù)據(jù)具有更高的信噪比(Martinez et al.，2013).隨著重力梯度全張量(FTG)儀器的問(wèn)世，F(xiàn)TG數(shù)據(jù)能夠提供更加豐富的多分量信息，每個(gè)梯度分量包含的信息量不盡相同，綜合利用各個(gè)張量信息有助于提高地質(zhì)解釋的準(zhǔn)確性.

隨著重力梯度測(cè)量的商業(yè)化和計(jì)算機(jī)技術(shù)的快速發(fā)展，重力梯度數(shù)據(jù)的三維反演技術(shù)成為國(guó)內(nèi)外學(xué)者研究的熱點(diǎn).重磁及其梯度三維反演問(wèn)題在數(shù)學(xué)上呈病態(tài).求解病態(tài)問(wèn)題可以得出多個(gè)解，這些解都滿(mǎn)足數(shù)據(jù)本身的要求，在地球物理反演中，被稱(chēng)為多解性問(wèn)題.引起多解性問(wèn)題主要有兩個(gè)原因：一方面，實(shí)際觀測(cè)數(shù)據(jù)是離散且有限的，并伴有來(lái)自?xún)x器和周?chē)h(huán)境的噪聲干擾，另一方面是場(chǎng)源等效性.通過(guò)增加野外觀測(cè)數(shù)據(jù)和提高觀測(cè)精度的辦法，可以在一定程度上降低反演的多解性，然而針對(duì)實(shí)際有限信息的觀測(cè)數(shù)據(jù)，更有效的方法是在反演中加入先驗(yàn)信息進(jìn)行約束反演(姚長(zhǎng)利等，2002).

多年來(lái)，約束條件的使用越來(lái)越受到重視，許多地球物理學(xué)者提出了各種不同的約束反演方法.Li和Oldenburg(1998)采用Tikhonov正則化方法對(duì)重力數(shù)據(jù)進(jìn)行三維反演，通過(guò)在反演目標(biāo)函數(shù)中引入粗糙度矩陣得到光滑反演模型，該方法也被稱(chēng)為光滑約束反演.Li(2001)將該方法應(yīng)用于重力梯度數(shù)據(jù)三維反演.Portniaguine和Zhdanov(1999)提出了聚焦反演方法并將其應(yīng)用于磁場(chǎng)數(shù)據(jù)、重力及重力梯度數(shù)據(jù)反演(Zhdanov et al.，2004)，反演結(jié)果表明該方法對(duì)具有尖銳邊界的地質(zhì)體有較好的反演效果.在國(guó)內(nèi)，姚長(zhǎng)利等(2007)針對(duì)三維物性反演中存在的多解性和內(nèi)存消耗巨大的問(wèn)題，提出了隨機(jī)子域反演技術(shù).該方法通過(guò)概率方式控制子域生成的分布，實(shí)現(xiàn)三維反演約束新機(jī)制.郭良輝等(2009)提出了基于重力梯度數(shù)據(jù)異常分離的三維相關(guān)成像方法.孟小紅等(2012)提出基于剩余異常相關(guān)成像的重磁物性反演方法，通過(guò)對(duì)物性模型的正演和實(shí)測(cè)結(jié)果的殘差進(jìn)行相關(guān)成像并迭代更新物性模型，實(shí)現(xiàn)了相關(guān)成像方法對(duì)物性參數(shù)的反演.

約束條件的使用在地球物理反演中至關(guān)重要，但約束條件的提取以及同反演計(jì)算的結(jié)合仍然存在很多困難，例如如何將一些地質(zhì)信息轉(zhuǎn)化成具體數(shù)學(xué)形式的約束條件還存在著困難，另外，有些約束條件還難以同反演計(jì)算過(guò)程融合，從而制約了約束條件作用的發(fā)揮(姚長(zhǎng)利等，2002；Li,2012).因此，發(fā)展新的容易與先驗(yàn)信息相結(jié)合的反演方法將具有重要實(shí)用價(jià)值(Sun and Li，2000).近年來(lái)，地質(zhì)統(tǒng)計(jì)學(xué)方法被應(yīng)用到地球物理反演中.地質(zhì)統(tǒng)計(jì)學(xué)方法在地球物理反演中應(yīng)用的好處，除了對(duì)區(qū)域地質(zhì)特征規(guī)律信息的應(yīng)用之外，還體現(xiàn)在能夠靈活地結(jié)合地質(zhì)學(xué)家對(duì)地質(zhì)情況的認(rèn)識(shí)上.Asli等(2000)建立了以待反演密度作為主變量，重力數(shù)據(jù)作為次級(jí)變量的兩變量協(xié)同克里金反演方程；Gloaguen等(2005,2007)將協(xié)同克里金反演方法應(yīng)用于井中雷達(dá)速度反演和層析成像領(lǐng)域；Shamsipour等(2010)采用協(xié)同克里金方法對(duì)重力數(shù)據(jù)進(jìn)行反演;Geng等(2014)采用協(xié)同克里金方法對(duì)重力梯度全張量數(shù)據(jù)進(jìn)行反演，反演結(jié)果表明重力梯度全張量數(shù)據(jù)聯(lián)合反演有助于提高反演結(jié)果的準(zhǔn)確性與分辨率，但反演結(jié)果在縱向上仍然存在一定程度的發(fā)散.

為了進(jìn)一步提高協(xié)同克里金反演方法結(jié)果的分辨率，提高該方法與先驗(yàn)信息有效融合的能力，本文在Geng等人(2014)工作基礎(chǔ)上，提出密度約束下的協(xié)同克里金反演方法.為檢驗(yàn)方法的有效性，研究密度信息的加入對(duì)反演結(jié)果的影響，對(duì)直立長(zhǎng)方體組合模型和傾斜巖脈組合模型的重力梯度全張量數(shù)據(jù)進(jìn)行了反演計(jì)算.最后將本文研究方法應(yīng)用于美國(guó)德克薩斯州的一個(gè)研究構(gòu)造上所獲得的數(shù)據(jù)，證明了本文方法的有效性和應(yīng)用前景.

2重力梯度反演的基本原理

2.1重力梯度張量與正演

在笛卡爾坐標(biāo)系下，對(duì)于地下任意形狀的地質(zhì)體已知其密度分布，其重力梯度張量，即重力場(chǎng)三分量沿3個(gè)坐標(biāo)系方向的梯度可表示為

(1)

由于在地球外部，引力位滿(mǎn)足拉普拉斯方程Txx+Tyy+Tzz=0，以及重力場(chǎng)的無(wú)旋性：Txy=Tyx，Txz=Tzx，Tyz=Tzy，因此，上式(1)中重力梯度張量的9個(gè)分量中只有5個(gè)是獨(dú)立的.將待反演的地下空間離散成一系列大小相等緊密排列的直立長(zhǎng)方體單元，并賦予每個(gè)單元固定的密度值，地表重力梯度張量各分量觀測(cè)值與地下密度分布可以表示為如下線(xiàn)性關(guān)：

(2)

其中m為網(wǎng)格單元的剩余密度，d為觀測(cè)點(diǎn)上的重力梯度張量分量的理論觀測(cè)值，Gα β(α,β=x,y,z)為該張量的正演核函數(shù)，它是一個(gè)離散單元相對(duì)觀測(cè)點(diǎn)位置所決定的常數(shù)矩陣(郭志宏等，2004).

2.2協(xié)同克里金反演

+ΛTCTα β,Tα βΛ,

(3)

其中， Λ是加權(quán)系數(shù)矩陣, Cρ ρ是密度協(xié)方差矩陣，CTα β,ρ是重力梯度數(shù)據(jù)和密度的互協(xié)方差矩陣, CTα β,Tα β是重力梯度數(shù)據(jù)的協(xié)方差矩陣(Aslietal., 2000).通過(guò)使主變量ρ估計(jì)方差最小，可以得到單分量反演的協(xié)同克里金方程：

(4)

根據(jù)式(4)，我們可以得到最優(yōu)加權(quán)系數(shù)矩陣Λ，然后利用最優(yōu)加權(quán)系數(shù)矩陣，我們可以得到密度估計(jì)值：

(5)

協(xié)同克里金方差向量可由下式得到：

(6)

估計(jì)方差位于式(6)的主對(duì)角線(xiàn)上，非對(duì)角元素是估計(jì)誤差的互協(xié)方差.由式(2)，我們知道密度和重力梯度數(shù)據(jù)間是線(xiàn)性關(guān)系，因此密度和重力梯度數(shù)據(jù)互協(xié)方差矩陣也是線(xiàn)性相關(guān)的(Myers, 1982)：

(7)

其中Q0是重力梯度觀測(cè)誤差協(xié)方差矩陣.觀測(cè)數(shù)據(jù)經(jīng)過(guò)一系列校正和處理后，可以認(rèn)為數(shù)據(jù)中的噪聲是不相關(guān)的，則矩陣Q0=σ2I.σ是觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)偏差，I為單位矩陣.

重力梯度數(shù)據(jù)和密度的互協(xié)方差與密度協(xié)方差之間存在如下的關(guān)系(Asli et al., 2000)：

(8)

不考慮重力梯度數(shù)據(jù)中的觀測(cè)誤差，即Q0=0，可以證明下式成立：

=Gα β(Gα βCρ ρ)T(Gα βCρ ρ(Gα β)T)-1Tα β=Tα β.

(9)

式(9)證明了采用協(xié)同克里金方法反演得到的模型的正演數(shù)據(jù)與觀測(cè)異常在不考慮噪聲影響的情況下是相等，從而保證了反演結(jié)果的正演數(shù)據(jù)與實(shí)際觀測(cè)數(shù)據(jù)能夠在噪聲水平內(nèi)得到較好的擬合.

如果已知先驗(yàn)信息，例如巖石的密度信息，用ρF表示，則建立密度數(shù)據(jù)約束下的重力梯度全張量數(shù)據(jù)反演的協(xié)同克里金方程如下式所示：

(10)

其中Λ,Γ,Ω,γ,H,Φ分別為分量Txz, Tyz, Tzz, Txx, Txy,ρF對(duì)應(yīng)的最優(yōu)加權(quán)矩陣.左側(cè)矩陣中對(duì)角元素為次級(jí)變量的協(xié)方差矩陣，非對(duì)角元素為不同次級(jí)變量之間的互協(xié)方差矩陣.對(duì)式(10)求解可以得到最優(yōu)加權(quán)矩陣，從而可以得到6個(gè)次級(jí)變量聯(lián)合反演的密度估計(jì)值：

ρ*=ΛTTxz+ΓTTyz+ΩTTzz+γTTxx+HTTxy+ΦTρF,

(11)

協(xié)同克里金估計(jì)方差可由下式計(jì)算得到：

(12)

由式(11)可知，采用協(xié)同克里金方法反演密度必須先計(jì)算出密度協(xié)方差矩陣Cρ ρ.要想得到Cρ ρ，需要先計(jì)算變差函數(shù)參數(shù)，即塊金值C0，x、y、z三方向的變程ax、ay、az，以及基臺(tái)值C+C0(C稱(chēng)為拱高).然后由計(jì)算得到的變差函數(shù)參數(shù)計(jì)算出密度協(xié)方差矩陣Cρρ(孫洪泉，1990).估計(jì)變差函數(shù)參數(shù)常用的方法有三種(Chasseriau and Chouteau，2003；Shamsipour et al.，2010)：(1)根據(jù)已有的地質(zhì)信息，建立初步估計(jì)模型，然后利用地質(zhì)統(tǒng)計(jì)學(xué)方法得到變差函數(shù)參數(shù)(附錄A)；(2)當(dāng)來(lái)自地下密度數(shù)據(jù)(例如測(cè)井密度)充足時(shí)，可以直接根據(jù)密度數(shù)據(jù)通過(guò)地質(zhì)統(tǒng)計(jì)學(xué)方法計(jì)算變差函數(shù)參數(shù)(附錄A)；(3)在先驗(yàn)地質(zhì)信息缺乏的情況下，可采用Alis等(2000)提出的V-V繪制方法估計(jì)變差函數(shù)參數(shù)(附錄B).

2.3積分靈敏度矩陣

在位場(chǎng)反演中核函數(shù)隨積分單元與觀測(cè)位置之間的距離增大而迅速衰減的性質(zhì)會(huì)導(dǎo)致反演結(jié)果的物性參數(shù)分布全部集中于地表，這種現(xiàn)象被稱(chēng)為反演的“趨膚效應(yīng)”.為了克服這種效應(yīng)，Li和Oldenburg(1996, 1998)在三維磁化率和密度反演中引入深度加權(quán)函數(shù)w(z)=(z+z0)-β.但該深度加權(quán)函數(shù)僅考慮深度方向的情況，未能考慮反演的橫向分辨率.此外，反演結(jié)果受參數(shù)β影響很大.因此本文在協(xié)同克里金方程中引入積分靈敏度矩陣.由于積分靈敏度矩陣能夠?qū),y,z三個(gè)方向加權(quán)，使得反演結(jié)果具有三個(gè)方向的分辨能力，從而使反演算法的橫向分辨能力得到提高.

觀測(cè)數(shù)據(jù)與模型變化的關(guān)系可表示為

(13)

其中，Gik為重力梯度數(shù)據(jù)的核函數(shù).數(shù)據(jù)靈敏度對(duì)模型mk的積分可以表示為

(14)

可以看到，數(shù)據(jù)靈敏度只與核函數(shù)有關(guān)，不涉及其他參數(shù)，從而避免了因參數(shù)選擇不同對(duì)反演結(jié)果的影響.在離散化的協(xié)同克里金反演方程中，積分靈敏度矩陣是一個(gè)對(duì)角矩陣：

W=diag(1/Sk).

(15)

向式(5)中引入積分靈敏度矩陣后，單分量協(xié)同克里金反演方程變?yōu)?/p>

(16)

3理論模型試驗(yàn)

為了檢驗(yàn)所建立的反演方程的正確性，并研究密度信息的加入對(duì)反演結(jié)果的影響，設(shè)計(jì)了兩個(gè)理論模型進(jìn)行反演試驗(yàn).理論模型試驗(yàn)的具體過(guò)程為：首先根據(jù)設(shè)計(jì)的理論模型采用快速正演算法(姚長(zhǎng)利等，2003)計(jì)算出重力梯度各個(gè)分量的正演數(shù)據(jù)，并向各個(gè)分量中加入5%的高斯白噪聲；然后對(duì)重力梯度張量的五個(gè)獨(dú)立分量進(jìn)行聯(lián)合反演，最后在已知的密度數(shù)據(jù)約束下對(duì)五個(gè)獨(dú)立分量進(jìn)行聯(lián)合反演.

3.1雙直立長(zhǎng)方體模型

首先設(shè)計(jì)由兩個(gè)直立長(zhǎng)方體組成的組合模型一.設(shè)置地下模型空間尺寸x，y，z為2400 m×2400 m×1500 m.將地下模型空間剖分為24×24×15個(gè)緊密排列的直立立方體單元，每個(gè)單元為100 m×100 m×100m.地面觀測(cè)網(wǎng)格為24×24的網(wǎng)格，觀測(cè)點(diǎn)位于每個(gè)網(wǎng)格中心.設(shè)地下模型空間中有兩個(gè)400 m×400 m×500 m的長(zhǎng)方體模型，剩余密度分別為-1 g·cm-3和1 g·cm-3.圖1a顯示了合成模型的密度分布三維圖.反演計(jì)算前，向該理論模型在地表產(chǎn)生的重力梯度張量的6個(gè)分量各加入最大幅值5%的高斯白噪聲.加入噪聲干擾后的Tzz分量如圖1b所示.

模型試驗(yàn)中，我們采用上述第一種方法計(jì)算變差函數(shù)參數(shù)值.采用高斯模型對(duì)估計(jì)模型(這里為真實(shí)模型)擬合結(jié)果如圖2(a—c)所示,得到變差函數(shù)參數(shù)：C0=2000 (kg·m-3)2，C=24000 (kg·m-3)2，ax=450 m,ay=450 m，az=450 m.

首先對(duì)張量的5個(gè)獨(dú)立分量Txy, Txz, Tyy, Tyz和Tzz進(jìn)行聯(lián)合協(xié)同克里金反演，反演結(jié)果在z=600 m橫向切片和y=1300 m縱向切片如圖3a和3b所示，其中黑色的方框表示真實(shí)模型的邊界.可以看到，反演結(jié)果較準(zhǔn)確地反映出真實(shí)異常體的空間位置和形態(tài)，但是剩余密度值的絕對(duì)值比真實(shí)值小，并且異常體底部密度值存在一定程度的發(fā)散.然后加入兩口測(cè)井中的密度數(shù)據(jù)，采用新建立的式(11)對(duì)五個(gè)獨(dú)立分量數(shù)據(jù)進(jìn)行聯(lián)合反演，其中兩口測(cè)井的水平中心坐標(biāo)分別為(0.8 km，1.3 km)和(1.8 km，1.3 km)，即測(cè)井恰好垂直穿過(guò)兩個(gè)異常體的中心，地表下方測(cè)井所穿過(guò)的網(wǎng)格單元的剩余密度為已知，模擬多井資料.圖3c和3d分別顯示了在密度約束下五個(gè)獨(dú)立分量聯(lián)合反演結(jié)果的橫向和縱向切片.可以看到與無(wú)密度約束反演結(jié)果相比，加入密度信息約束反演得到的異常體邊界更加清晰，尤其是縱向分辨率得到顯著提高；此外，還可以看到測(cè)井穿過(guò)的單元附近的剖分單元也能夠比較準(zhǔn)確反演.加入密度數(shù)據(jù)前后反演結(jié)果與真實(shí)模型的均方根誤差分別為0.10 g·cm-3和0.07 g·cm-3，從而從定量角度說(shuō)明了加入密度數(shù)據(jù)進(jìn)行約束反演能夠提高反演結(jié)果的精度.

圖1　(a) 組合長(zhǎng)方體理論模型; (b) 理論模型Tzz分量正演結(jié)果(含5%的高斯隨機(jī)噪聲)Fig.1　(a) The synthesized prism model; (b) calculated Tzz component with 5% uncorrelated Gaussian noise of maximum of the datum magnitude

圖2　高斯模型擬合變差函數(shù)(a) x方向擬合結(jié)果； (b) y方向擬合結(jié)果； (c) z方向擬合結(jié)果.Fig.2　Fitting variogram function using Gaussian model(a) x-direction variogram; (b) y-direction variogram; (c) z-direction variogram.

圖3　無(wú)密度約束下張量的五個(gè)獨(dú)立分量聯(lián)合反演結(jié)果：(a) z=600 m橫向切片； (b) y=1200 m縱向切片.密度約束下張量的五個(gè)獨(dú)立分量聯(lián)合反演結(jié)果圖： (c) z=600 m橫向切片； (d) y=1200 m縱向切片F(xiàn)ig.3　The recovered model with five independent componentswithout constraints： (a) In the z=600 m depth section; (b) In the y=1200 m cross section. The recovered model with five independent components with constraints： (c) In the z=600 m depth section; (d) In the y=1200 m cross section

圖4　(a) 組合傾斜脈狀體理論模型; (b)理論模型Tzz分量正演結(jié)果(含5%的高斯隨機(jī)噪聲)Fig.4　(a) The synthesized dike model; (b) Calculated Tzz component with 5% uncorrelated Gaussian noise of the datum magnitude

圖5　高斯模型擬合變差函數(shù)(a) y方向擬合結(jié)果； (b) 平行于脈狀體傾向方向擬合結(jié)果； (c) 垂直于脈狀體傾向方向擬合結(jié)果.Fig.5　Fitting variogram function using Gaussian model(a) y-direction variogram; (b) Variogram along the inclination of dikes; (c) Variogram perpendicular to the inclination of dikes.

圖6　無(wú)密度約束下張量的五個(gè)獨(dú)立分量聯(lián)合反演結(jié)果： (a) z=600 m橫向切片; (b) y=1200 m縱向切片.密度約束下張量的五個(gè)獨(dú)立分量聯(lián)合反演結(jié)果： (c) z=600 m橫向切片; (d) y=1200 m縱向切片F(xiàn)ig.6　The recovered model with five independent componentswithout constraints： (a) The z=600 m depth section; (b) The y=1200 m cross section. The recovered model with five independent components with constraints: (c) The z=600 m depth section; (d) The y=1200 m cross section

圖7　重力梯度全張量異常(a) Txx; (b) Txy; (c) Tzx; (d) Tyy; (e) Tyz; (f) Tzz.Fig.7　Gravity gradient tensor data

圖8　重力梯度張量的五個(gè)獨(dú)立分量聯(lián)合反演結(jié)果： (a) 剩余密度分布三維圖; (b) z=4500橫向切片; (c) 巖丘走向方向縱向切片; (d) 垂直于走向方向縱向切片. 密度約束下五個(gè)獨(dú)立分量聯(lián)合反演結(jié)果： (e) 剩余密度分布三維圖; (f) z=4500 m橫向切片; (g) 巖丘走向方向縱向切片; (h) 垂直于走向方向縱向切片F(xiàn)ig.8　The recovered model with five independent components： (a) Recovered 3-D model; (b) The depth slice of 4500 m; (c) Cross section along the strike direction of salt dome; (d) Cross section perpendicular to the strike direction of salt dome. The recovered model with five independent components with constraints: (e) Recovered 3-D model; (f) The depth an>slice of 4500 m; (g) Cross section along the strike direction of salt dome; (h) Cross section perpendicular to the strike direction of salt dome

3.2脈狀體組合模型

下面設(shè)計(jì)一個(gè)較為復(fù)雜的傾斜脈狀體組合模型，如圖4a所示，兩個(gè)脈狀體的傾角都為45°，剩余密度分別為0.8 g·cm-3和1 g·cm-3.我們知道，位場(chǎng)數(shù)據(jù)的縱向分辨率非常差，縱向上存在疊加的異常體很難被區(qū)分.這一組合模型用來(lái)研究密度約束在反演復(fù)雜模型中的作用.設(shè)置地下模型空間尺寸x,y,z為2400 m× 2400 m ×1500 m.將地下模型空間剖分為24×24×15個(gè)緊密排列的直立立方體單元，每個(gè)單元為100 m×100 m×100 m.地面觀測(cè)網(wǎng)格為24×24的網(wǎng)格，觀測(cè)點(diǎn)位于每個(gè)網(wǎng)格中心.反演計(jì)算前，向該理論模型在地表產(chǎn)生的重力梯度張量的六個(gè)分量各加入5%的高斯白噪聲，帶有噪聲干擾的Tzz分量如圖4b所示.

采用高斯模型對(duì)估計(jì)模型(這里為真實(shí)模型)進(jìn)行擬合.擬合結(jié)果如圖5(a—c)所示，得到變差函數(shù)參數(shù)：C0=2000 (kg·m-3)2，C=24000 (kg·m-3)2，ay=250 m，a45°=1000 m，a135°=380 m，其中a45°的方向?yàn)槠叫杏诿}狀體傾向的方向，a135°的方向垂直于y軸和脈狀體傾向方向所構(gòu)成的平面.首先對(duì)5個(gè)獨(dú)立分量Txy，Txz，Tyy，Tyz和Tzz進(jìn)行聯(lián)合協(xié)同克里金反演，反演結(jié)果在z=600 m橫向切片和y=1300 m縱向切片如圖6a和6b所示，可以看到縱向上存在疊加的兩個(gè)傾斜脈狀體都得到了比較準(zhǔn)確的恢復(fù)：兩個(gè)脈狀體的位置、傾角都與真實(shí)模型基本一致，但是脈狀體底部密度值發(fā)散比較嚴(yán)重.然后加入一口測(cè)井中的密度數(shù)據(jù)作為約束條件再對(duì)五個(gè)獨(dú)立分量數(shù)據(jù)進(jìn)行聯(lián)合反演，測(cè)井的水平中心坐標(biāo)為(1.3 km，1.3 km)，地表下方測(cè)井所穿過(guò)的網(wǎng)格剩余密度已知，模擬單口井資料.圖6c和6d分別顯示了五個(gè)獨(dú)立分量在密度約束下聯(lián)合反演結(jié)果的橫向和縱向切片圖，可以看到與無(wú)密度約束反演結(jié)果相比，加入密度信息約束后，反演得到的兩個(gè)傾斜脈狀體的邊緣發(fā)散現(xiàn)象得到明顯改善，脈狀體邊緣更加清晰，并且兩個(gè)脈狀體的剩余密度與各自的真實(shí)值也非常接近.加入密度數(shù)據(jù)前后反演結(jié)果與真實(shí)模型的均方根誤差分別為0.09 g·cm-3和0.07 g·cm-3，再次證明了密度數(shù)據(jù)的加入對(duì)提高重力梯度全張量數(shù)據(jù)反演精度的能力.

4實(shí)測(cè)數(shù)據(jù)應(yīng)用

將協(xié)同克里金反演方法應(yīng)用于美國(guó)德克薩斯州的一個(gè)研究構(gòu)造上所獲得的重力數(shù)據(jù)(Nettleton，1976)

的解釋中.利用FFT變換(Michus and Hinojosa，2001)計(jì)算得到異常的全張量數(shù)據(jù).研究區(qū)內(nèi)為沉積巖相地下發(fā)育鹽丘構(gòu)造.石油和天然氣的沉積及儲(chǔ)層與鹽丘有密切關(guān)系，鹽丘區(qū)附近常常擁有豐富的油氣儲(chǔ)量，因此對(duì)鹽丘地區(qū)的研究具有重要價(jià)值(Zhou，2006；Ennen and Hall，2011；梁杰等，2010).由于鹽丘下部構(gòu)造往往非常復(fù)雜，地震波在鹽丘下部迅速衰減，地震照明度低，其地震成像面臨著非常大的挑戰(zhàn)(Liu et al.,2011).在鹽丘地區(qū)，對(duì)重力或重力梯度全張量數(shù)據(jù)進(jìn)行三維反演，得到的密度模型可為該地區(qū)的研究工作提供依據(jù)；此外，得到的密度模型可以通過(guò)適當(dāng)?shù)拿芏?速度函數(shù)轉(zhuǎn)化成速度模型，從而實(shí)現(xiàn)地震-重力資料的綜合解釋?zhuān)岣呓忉尳Y(jié)果的可靠性.

重力梯度反演選擇重力梯度異常低值區(qū)域13 km×13 km，并將數(shù)據(jù)網(wǎng)格化，網(wǎng)格間距是500 m，網(wǎng)格化后的重力梯度張量的各個(gè)分量等值線(xiàn)如圖7所示.將反演目標(biāo)區(qū)域地下空間劃分為26×26×20個(gè)緊密排列的直立正方體單元，每個(gè)單元的尺寸為500 m×500 m×500 m.采用V-V繪制方法得到變差函數(shù)參數(shù)值:C0=2000 (kg·m-3)2，C=26000 (kg·m-3)2，a-40°=5000 m(巖丘延伸方向)，a50°=4000 m(垂直巖丘延伸方向)，az=3000 m.

首先對(duì)五個(gè)獨(dú)立分量Txx，Txy，Txz，Tyz，Tzz進(jìn)行聯(lián)合協(xié)同克里金反演，反演結(jié)果剩余密度三維分布如8a所示，反演結(jié)果在z=4500 m的橫向切片如圖8b所示，巖丘延伸方向的縱向切片(圖8b上虛線(xiàn)l1所示的位置)和垂直巖丘延伸方向的縱向切片(圖8b虛線(xiàn)l2所示的位置)分別如圖8c和8d所示.為了檢驗(yàn)密度數(shù)據(jù)約束反演在提高反演結(jié)果分辨率方面的能力，設(shè)置了五口測(cè)井，已有地質(zhì)調(diào)查和鉆井?dāng)?shù)據(jù)顯示(Said,1962)該巖丘的剩余密度在-0.3 g·cm-3左右，因此令這五口測(cè)井穿越過(guò)的巖石的剩余密度為-0.3 g·cm-3.在五口井中密度數(shù)據(jù)約束下，對(duì)重力梯度全張量數(shù)據(jù)進(jìn)行協(xié)同克里金反演，反演結(jié)果剩余密度三維分布如圖8e所示，相對(duì)應(yīng)的切片結(jié)果如圖8(f、g、h)所示(其中白色的線(xiàn)和點(diǎn)表示測(cè)井所在位置和穿越的單元格).對(duì)比圖8(a—d)和8(e—h)可以明顯看到，在密度數(shù)據(jù)的約束下，反演結(jié)果的分辨率得到明顯提高，測(cè)井穿過(guò)的單元附近的剖分單元剩余密度也非常接近-0.3 g·cm-3.從鹽丘區(qū)的反演結(jié)果看，低密度值分布約在地下3～6.5 km的深度；巖丘中心埋深約為4.5 km，這一結(jié)果與以往研究結(jié)果(見(jiàn)表1)相一致.本文采用的5組密度數(shù)據(jù)是根據(jù)已有調(diào)查和研究結(jié)果所模擬的數(shù)據(jù)，在實(shí)際數(shù)據(jù)應(yīng)用中，在缺少井?dāng)?shù)據(jù)的情況下，可以將地質(zhì)調(diào)查和/或其他物探方法得到的密度信息作為約束進(jìn)行聯(lián)合反演，同樣可以起到降低反演多解性，提高反演結(jié)果的分辨率的作用.

表1　密度約束下協(xié)同克里金反演結(jié)果與先前計(jì)算巖丘中心埋深結(jié)果比較

5結(jié)論

本文建立了密度約束下的重力梯度全張量數(shù)據(jù)聯(lián)合協(xié)同克里金反演方程，使協(xié)同克里金反演方法不僅可以將地質(zhì)學(xué)家對(duì)地質(zhì)情況的認(rèn)識(shí)，例如地質(zhì)體走向、傾角、空間尺度等結(jié)合到反演方程，還能夠?qū)?lái)自于測(cè)井或其他方法獲得的密度數(shù)據(jù)結(jié)合到反演方程中，從而大大降低反演的多解性.理論模型試驗(yàn)表明加入密度數(shù)據(jù)進(jìn)行約束反演能夠顯著提高反演結(jié)果的分辨率與精度.最后將其應(yīng)用于美國(guó)德克薩斯州鹽丘的重力梯度數(shù)據(jù)的解釋?zhuān)@得結(jié)果與前人解釋結(jié)果相一致，從而證明了該方法的實(shí)用性.

致謝特別感謝提供實(shí)際資料的同仁以及為評(píng)審本文所付出努力的專(zhuān)家們．

附錄A變差函數(shù)的估計(jì)與擬合

在地質(zhì)統(tǒng)計(jì)學(xué)中，變差函數(shù)是最基本也是最重要的模擬工具，它用于描述數(shù)據(jù)值的空間互相關(guān)性，將數(shù)據(jù)點(diǎn)的相關(guān)性隨空間距離的增大而減小的現(xiàn)象通過(guò)數(shù)學(xué)函數(shù)模擬出來(lái).變差函數(shù)可以用四個(gè)參數(shù)來(lái)描述，即變差函數(shù)類(lèi)型、變程a、塊金效應(yīng)C0和基臺(tái)值C+C0.

(1)變差函數(shù)的計(jì)算公式與估算

變差函數(shù)的定義是：區(qū)域化變量Z(x)和Z(x+h)兩點(diǎn)之差的方差之半.Z(x)的變差函數(shù)數(shù)學(xué)定義如下(孫洪泉，1990)：

(A1)

其中Var(·)表示方差，h為滯后距.如果有了區(qū)域變量Z(x)的一部分采樣樣本，就可以估算該區(qū)域變化量Z(x)的變差函數(shù)，具體計(jì)算公式如下：

(A2)

其中i為樣本序號(hào).當(dāng)來(lái)自地下(例如通過(guò)測(cè)井方法得到)密度數(shù)據(jù)量充足時(shí)，變差函數(shù)可以直接利用已知的密度數(shù)據(jù)根據(jù)式(A2)計(jì)算得到；或解釋人員根據(jù)已有的地質(zhì)信息，將地下解釋空間剖分成緊密排列的直立棱柱體單元，建立初步解釋模型，然后再根據(jù)式(A2)計(jì)算變差函數(shù).

(2)變差函數(shù)的擬合

下面我們給出根據(jù)初步解釋模型估算變差函數(shù)的示例.

設(shè)置地下模型空間尺寸x，y，z為2100 m× 2100 m ×1500 m.將地下模型空間剖分為21×21×15個(gè)緊密排列的直立立方體單元，每個(gè)單元尺寸為100 m×100 m×100 m.設(shè)在地下模型空間中有一個(gè)300 m×500 m×700 m的長(zhǎng)方體模型，剩余密度為1 g·cm-3，如圖A1所示.

圖A1　直立長(zhǎng)方體模型Fig.A1　Prism model

圖A2　高斯模型擬合變差函數(shù)Fig.A2　Fitting variogram function using Gaussian model

由于這一地質(zhì)模型屬性具有各向異性特征，應(yīng)當(dāng)計(jì)算它在不同方向的變差函數(shù)值.滯后距h是指向某方向移動(dòng)的距離，計(jì)算時(shí)，在指定的方向上，對(duì)指定的h，搜索所有相距h的點(diǎn)對(duì)[Z(xi), [Z(xi+h)]，并統(tǒng)計(jì)點(diǎn)對(duì)個(gè)數(shù)N，根據(jù)式(A2)計(jì)算出滯后距h對(duì)應(yīng)的變差值.然后再計(jì)算出滯后距2h、3h、4h…對(duì)應(yīng)的變差值.圖A2(a—c)中的小圓圈分別給出了x，y，z三個(gè)方向的變差的計(jì)算結(jié)果.在統(tǒng)計(jì)模擬中，變差函數(shù)采用的是一個(gè)解析函數(shù)，因此，計(jì)算出不同滯后距h的變差值V(h)后，還需要擬合出一個(gè)解析函數(shù).變差函數(shù)的擬合有很多種方法，比如在指定出函數(shù)類(lèi)型后，可以使用加權(quán)最小二乘法、線(xiàn)性規(guī)劃法、遺傳算法等.其中遺傳算法能夠擬合任何指定的函數(shù)類(lèi)型，因此，本文采用遺傳算法進(jìn)行擬合.

常用的解析函數(shù)類(lèi)型有三種，分別是球狀模型、高斯模型和指數(shù)模型.選用哪種函數(shù)擬合，一方面取決于樣本數(shù)據(jù)，另一方面，也取決于對(duì)實(shí)際規(guī)律的掌握程度.這里我們采用高斯模型，它的形式是

(A3)

附錄BV-V繪制方法

V-V繪制方法實(shí)現(xiàn)步驟：

(1) 選擇一個(gè)變差模型，包括解析函數(shù)類(lèi)型和三個(gè)參數(shù)：基臺(tái)值、變程和塊金值.由此，計(jì)算出密度協(xié)方差矩陣Cρ ρ.

(3) 對(duì)理論協(xié)方差從小到大進(jìn)行排序并分組，用同樣的排列順序?qū)?shí)驗(yàn)協(xié)方差進(jìn)行排序并分組.

(4) 計(jì)算分完組后的理論協(xié)方差和實(shí)驗(yàn)協(xié)方差的平均絕對(duì)誤差.

(5) 修改變差模型參數(shù)使理論協(xié)方差和實(shí)驗(yàn)協(xié)方差的平均絕對(duì)誤差最小，從而使二者之間具有較高的相關(guān)性.擬合過(guò)程中，結(jié)合已有的地質(zhì)信息修改變差模型參數(shù)，有利于得到好的擬合結(jié)果.

References

Abdelrahman E M, Bayoumi A I, Elaraby H M.Short note: A least-squares minimization approach to invert gravity data.Geophysics, 1991, 56: 115-118.

Asli M, Marcotte D, Chouteau M. 2000. Direct inversion of gravitydata by cokriging.∥Kleingeld W,Krige D eds.Geostats 2000, Proceedings of the 6th International Geostatistics Congress, 64-73.ChasseriauP, Chouteau M. 2003. 3D gravity inversion using a modelof parameter covariance.JournalofAppliedGeophysics, 52(1): 59-74.Chilès J P, Delfiner P. Geostatistics: Modeling Spatial Uncertainty, Second Edition, Hoboken, NJ, Wiley, 2012.

Ennen C, Hall S. 2011. Structural mapping of the Vinton salt dome, Louisiana, using gravity gradiometry data.∥ 2011SEGAnnual Meeting Technical Program. Society of Exploration Geophysicists.Expanded Abstracts, 830-835.

Essa K S.Gravity data interpretation using the s-curves method.J.Geophys.Eng., 2007, 4(2): 204-213.

Geng M X, Huang D N, Yang Q J, et al. 2014. 3D inversion of airborne gravity-gradiometry data using cokriging.Geophysics, 79(4): G37-G47.

Gloaguen E, Marcotte D, Chouteau M, et al. 2005. Borehole radar velocity inversion using cokriging and cosimulation.JournalofAppliedGeophysics, 57(4): 242-259.

Gloaguen E, Marcotte D, Giroux B, et al. 2007. Stochastic borehole radar velocity and attenuation tomographies using cokriging and cosimulation.JournalofAppliedGeophysics, 62(2): 141-157.

Guo L H, Meng X H, Shi L, et al. 2009. 3-D correlation imaging for gravity and gravity gradiometry data.ChineseJ.Geophys. (in Chinese), 52(4): 1098-1106, doi:10.3969/j.issn.0001-5733.2009.04.027.

Guo Z H, Guan Z N, Xiong S Q. 2004. Cuboid ΔT and its gradient forward theoretical expressions without analytic odd points.ChineseJ.Geophys.(in Chinese), 47(6): 1131-1138, doi:10.3321/j.issn:0001-5733.2004.06.029.

Li Y G, Oldenburg D W. 1998. 3-D inversion of gravity data.Geophysics, 63(1): 109-119.

Li Y G. 2001. 3-D inversion of gravity gradiometer data.∥ 2001 SEG Annual Meeting Technical Program.Society of Exploration Geophysicists.Expanded Abstracts,1470-1473. Li Y G. 2012. Recent advances in 3D generalized inversion of potential-field data. ∥ 2001 SEG Annual MeetingTechnical Program.Society of Exploration Geophysicists.Expanded Abstracts, 1-7.

Liang J, Gong J M, Cheng H Y. 2010. Control of salt rock distribution on oil and gas pooling in the Gulf of Mexico.MarineGeologyLetters(in Chinese), 26(1): 25-30.

Liu YK, Chang X, Jin D G, et al. 2011. Reverse time migration of multiples for subsalt imaging.Geophysics, 76(5): WB209-WB216.

Ma G Q, Du X J, Li L L. 2012. Comparison of the tensor local wavenumber method with the conventional local wavenumber method for interpretation of total tensor data of potential fields.ChineseJ.Geophys. (in Chinese), 55(7): 2450-2461, doi: 10.6038/j.issn.0001-5733.2012.07.029.

Martinez C, Li Y G, Krahenbuhl R,et al. 2013.3D inversion of airborne gravity gradiometry data in mineral exploration: A case study in the Quadrilátero Ferrífero, Brazil.Geophysics, 78(1): B1-B11.

Meng X H, Liu G F, Chen Z X, et al. 2012. 3-D gravity and magnetic inversion for physical properties based on residual anomaly correlation.ChineseJ.Geophys. (in Chinese), 55(1): 304-309, doi: 10.6038/j.issn.0001-5733.2012.01.030.

Mickus K L, Hinojosa J H. 2001. The complete gravity gradient tensor derived from the vertical component of gravity: A Fourier transform technique.JournalofAppliedGeophysics, 46(3): 159-174.

Mohan N L, Anandababu L, Roa S.Gravity interpretation using Mellin transform.Geophysics, 1986, 52(1): 114-122.

Myers D E. 1982.Matrix formulation of co-kriging.MathematicalGeology, 14(3): 249-257.

Nettleton L L. 1976. Gravity and Magnetic in Oil Prospecting. New York: McGraw-Hill Book Co.

Qruc B.Depth estimation of simple causative sources from gravity gradient tensor invariants and vertical component.PureandAppliedGeophysics, 2010, 167(10):1259-1272.

Portniaguine O, Zhdanov M S. 1999.Focusing geophysical inversion images.Geophysics, 64(3): 874-887.

Said R. 1962.The Geology of Egypt.Amsterdam: Elsevier.

Salem A, Ravatz D, Mushayandebvu M F, et al.Linearized least-squares method for interpretation of potential-field data from sources of simple geometry.Geophysics, 2004, 69(3): 783-788.

Shamsipour P, Marcotte D, Chouteau M, et al. 2010. 3D stochastic inversion of gravity data using cokriging and cosimulation.Geophysics, 75(1):I1-I10.

Shaw R K, Agarwal B N P.A generalized concept of resultant gradient to interpret potential field maps.GeophysicalProspecting, 1997, 45(3): 513-520.Sun H Q. 1990. Geological Statistics and Application. Beijing: China University of Mining and Technology Press.

Sun J J, Li Y G. 2010.Inversion of surface and borehole gravity with thresholding and density constraints.∥ 2010 SEG Annual Meeting Technical Program.Society of Exploration Geophysicists.Expanded Abstracts, 1798-1803.Yao C L, Hao T Y. Guan Z N. 2002.Restrictions in gravity and magnetic inversions and technical strategy of 3D properties inversion.Geophysical&GeochemicalExploration(in Chinese), 26(4): 253-257.Yao C L, Hao T Y, Guan Z N. 2003.High-speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms.ChineseJ.Geophys. (in Chinese), 46(2): 252-258, doi:10.3321/j.issn:0001-5733.2003.02.020.

Yao C L, Zheng Y M, Zhang Y W. 2007. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces.ChineseJ.Geophys. (in Chinese), 50(5): 1576-1583, doi:10.3321/j.issn:0001-5733.2007.05.035.

Zeng H L. 1999. Present state and revival of gravity gradiometry.Geophysical&GeochemicalExploration(in Chinese), 23(1): 1-6.

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

Zhou H W. 2006. First-break vertical seismic profiling tomography for Vinton Salt Dome.Geophysics, 71(3): U29-U36.

附中文參考文獻(xiàn)

郭良輝, 孟曉紅, 石磊等. 2009. 重力和重力梯度數(shù)據(jù)三維相關(guān)成像. 地球物理學(xué)報(bào), 52(4): 1098-1106, doi:10.3969/j.issn.0001-5733.2009.04.027.

郭志宏, 管志寧, 熊盛青. 2004. 長(zhǎng)方體ΔT場(chǎng)及其梯度場(chǎng)無(wú)解析奇點(diǎn)理論表達(dá)式. 地球物理學(xué)報(bào), 47(6): 1131-1138, doi:10.3321/j.issn:0001-5733.2004.06.029.

梁杰, 龔建明, 成海燕. 2010. 墨西哥灣鹽巖分布對(duì)油氣成藏的控制作用. 海洋地質(zhì)動(dòng)態(tài), 26(1): 25-30.

馬國(guó)慶, 杜曉娟, 李麗麗. 2012. 解釋位場(chǎng)全張量數(shù)據(jù)的張量局部波數(shù)法及其與常規(guī)局部波數(shù)法的比較. 地球物理學(xué)報(bào), 55(7): 2450-2461, doi: 10.6038/j.issn.0001-5733.2012.07.029.

孟小紅, 劉國(guó)峰, 陳召曦等. 2012. 基于剩余異常相關(guān)成像的重磁物性反演方法. 地球物理學(xué)報(bào), 55(1): 304-309, doi: 10.6038/j.issn.0001-5733.2012.01.030.

孫洪泉. 1990. 地質(zhì)統(tǒng)計(jì)學(xué)及其應(yīng)用.北京: 中國(guó)礦業(yè)大學(xué)出版社.

姚長(zhǎng)利, 郝天珧, 管志寧. 2002. 重磁反演約束條件及三維物性反演技術(shù)策略. 物探與化探, 26(4): 253-257.

姚長(zhǎng)利, 郝天珧, 管志寧等. 2003. 重磁遺傳算法三維反演中高速計(jì)算及有效存儲(chǔ)方法技術(shù). 地球物理學(xué)報(bào), 46(2): 252-258, doi:10.3321/j.issn:0001-5733.2003.02.020.

姚長(zhǎng)利, 鄭滿(mǎn)元, 張聿文. 2007. 重磁異常三維物性反演隨機(jī)子域法方法技術(shù). 地球物理學(xué)報(bào), 50(5): 1576-1583, doi:10.3321/j.issn:0001-5733.2007.05.035.

曾華霖. 1999. 重力梯度測(cè)量的現(xiàn)狀及復(fù)興. 物探與化探, 23(1): 1-6.

(本文編輯胡素芳)

Three-dimensional constrained inversion of full tensor gradiometer data based on cokriging method

GENG Mei-Xia1,2, HUANG Da-Nian2, YU Ping2*, YANG Qing-Jie3

1InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China2Geo-ExplorationScienceandTechnologyInstitute,JilinUniversity,Changchun130026,China3ChineseAcademyofSciences,InstituteofGeodesyandGeophysics,Wuhan430077,China

AbstractAs full tensor gradiometer (FTG) measurement techniques has increasingly become mature, three-dimensional inversion of FTG data is focused on more and more. Like inversion of gravity, inversion of FTG data confronts with multi-solution problem. In this paper, cokriging inversion equation is framed with several components with known densities as constrains. The method proposed can easily conclude the prior information, for example the known densities, into the inversion equation, so that non-uniqueness of inversion can be reduced. This method is tested on sythetic models and the inverted results show that the resolution, especially the vertical resolution of the resulted model can be significantly improved with the densities as constrains. The application of actual data obtained from a salt dome located in Texas proved the validity of the proposed method.

KeywordsFull tensor gradiometer； 3-D constrained inversion； Cokriging method； Geostatistics

基金項(xiàng)目國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)課題(2014AA06A613),中國(guó)博士后科學(xué)基金資助項(xiàng)目(2015M580680),中央高?；究蒲袠I(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金聯(lián)合資助.

作者簡(jiǎn)介耿美霞，女，講師，1986年生，河北石家莊人，主要從事重磁及其梯度數(shù)據(jù)反演研究. E-mail: gengmeixia@126.com E-mail: yuping@jlu.edu.cn

*通訊作者于平，女，博士，教授，1978年生，遼寧營(yíng)口人，主要從事綜合地球物理數(shù)據(jù)處理-解釋-建模及軟件集成一體化研究.

doi:10.6038/cjg20160528 中圖分類(lèi)號(hào)P631

收稿日期2015-09-02，2015-12-09收修定稿

耿美霞， 黃大年， 于平等. 2016. 基于協(xié)同克里金方法的重力梯度全張量三維約束反演.地球物理學(xué)報(bào)，59(5):1849-1860,doi:10.6038/cjg20160528.

Geng M X, Huang D N, Yu P, et al. 2016. Three-dimensional constrained inversion of full tensor gradiometer data based on cokriging method.ChineseJ.Geophys. (in Chinese),59(5):1849-1860,doi:10.6038/cjg20160528.