重力和重力梯度數(shù)據(jù)聯(lián)合聚焦反演方法

秦朋波， 黃大年

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春　130000

?

重力和重力梯度數(shù)據(jù)聯(lián)合聚焦反演方法

秦朋波， 黃大年

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春130000

摘要重力數(shù)據(jù)包含較多的低頻信息，重力梯度數(shù)據(jù)包含較多的高頻信息，將重力數(shù)據(jù)和重力梯度數(shù)據(jù)進(jìn)行聯(lián)合反演得到的結(jié)果更加可信.本文基于聚焦反演方法，實(shí)現(xiàn)了這一過(guò)程.因?yàn)槁?lián)合反演中分量種類增加，所以計(jì)算靈敏度矩陣所需要的時(shí)間增加，為此，本文提出了一種快速計(jì)算靈敏度矩陣的方法.因?yàn)槁?lián)合反演對(duì)內(nèi)存的要求增大，本文選擇有限內(nèi)存BFGS擬牛頓法求解反演問(wèn)題.本文通過(guò)再加權(quán)的方法實(shí)現(xiàn)深度加權(quán).文中利用單一分量的反演結(jié)果來(lái)預(yù)測(cè)異常體的埋深信息，隨后將埋深信息結(jié)合到深度加權(quán)函數(shù)中，將其用于多分量組合反演計(jì)算.給出了模型試驗(yàn)，發(fā)現(xiàn)預(yù)測(cè)得到的異常體的埋深信息與其實(shí)際埋深存在偏差，但是將這一信息應(yīng)用到反演計(jì)算，能夠得到與真實(shí)模型一致的結(jié)果.之后，本文通過(guò)模型試驗(yàn)來(lái)探究重力和重力梯度聯(lián)合反演的優(yōu)勢(shì)，發(fā)現(xiàn)將重力和重力梯度數(shù)據(jù)聯(lián)合，能夠識(shí)別出額外的噪聲，反演得到的模型更加合理.但是，對(duì)于不同分量組合得到的反演結(jié)果是相近的，反演模型的提高很小.最后，將聯(lián)合反演方法應(yīng)用到美國(guó)路易斯安那州Vinton巖丘的實(shí)際數(shù)據(jù)中，結(jié)果顯示，將重力和重力梯度數(shù)據(jù)聯(lián)合反演，反演模型得到了提高，反演得到的結(jié)果與地質(zhì)資料吻合.

關(guān)鍵詞重力和重力梯度數(shù)據(jù)正演； 重力和重力梯度聯(lián)合反演； 有限內(nèi)存BFGS擬牛頓法； 深度加權(quán)函數(shù)； 最小梯度支撐函數(shù)

1引言

重力數(shù)據(jù)反演和重力梯度數(shù)據(jù)反演方法是位場(chǎng)反演方法的重要組成部分.在早期，因?yàn)橹亓?shù)據(jù)獲取和處理更加容易，對(duì)計(jì)算機(jī)的要求低，人們提出了很多用于重力數(shù)據(jù)反演的方法(Last and Kubik，1983；Li and Oldenburg，1998；Li and Chouteau，1998；Portniaguine and Zhdanov，1999；Chasseriau and Chouteau，2003；Shamsipour and Marcotte，2010；劉銀萍等，2013).后來(lái)，隨著FTG系統(tǒng)的出現(xiàn)，梯度數(shù)據(jù)的精度得到了很大的提高，人們將很多反演方法延伸到了梯度數(shù)據(jù)的反演(Li，2001；Zhdanov et al.，2004；Martinez et al.，2013；Geng et al.，2014).在實(shí)際測(cè)量中，重力數(shù)據(jù)是通過(guò)測(cè)量位場(chǎng)的垂直分量獲得的，重力梯度數(shù)據(jù)則是通過(guò)測(cè)量重力場(chǎng)在三個(gè)方向的變化獲得，在頻率域?qū)煞N數(shù)據(jù)進(jìn)行比較，可以發(fā)現(xiàn)，重力數(shù)據(jù)包含較多的低頻信息，重力梯度數(shù)據(jù)則包含較多的高頻信息.基于這一點(diǎn)，將兩種數(shù)據(jù)結(jié)合進(jìn)行聯(lián)合反演得到的結(jié)果更加可信.但是，目前能夠?qū)⒅亓椭亓μ荻葦?shù)據(jù)聯(lián)合反演的方法還比較少.Zhdanov等(2004)通過(guò)重力曲率將重力數(shù)據(jù)和梯度數(shù)據(jù)結(jié)合進(jìn)行了反演計(jì)算.Wu等(2013)將目標(biāo)體看作是質(zhì)點(diǎn)，對(duì)其重力和重力梯度公式進(jìn)行轉(zhuǎn)換，實(shí)現(xiàn)了重力和重力梯度數(shù)據(jù)的聯(lián)合反演.Capriotti等(2014)針對(duì)重力和重力梯度數(shù)據(jù)核函數(shù)隨深度衰減速率不同的情況，通過(guò)靈敏度矩陣平衡兩種衰減速率，實(shí)現(xiàn)了重力和重力梯度數(shù)據(jù)的聯(lián)合反演.本文基于聚焦反演方法，提出一種能夠?qū)崿F(xiàn)重力和重力梯度數(shù)據(jù)聯(lián)合反演的方法.

進(jìn)行聯(lián)合反演計(jì)算之前，首先需要計(jì)算各分量的靈敏度矩陣.一般將地下規(guī)則劃分為一定數(shù)量的小長(zhǎng)方體，分量種類的增加和小長(zhǎng)方體數(shù)量的增加都會(huì)使靈敏度矩陣的計(jì)算時(shí)間增加.同時(shí)，存儲(chǔ)矩陣所需要的內(nèi)存也增加.Li等(2003)對(duì)靈敏度矩陣進(jìn)行小波變換，對(duì)于變換后的矩陣，通過(guò)設(shè)置閾值，將低于閾值的部分設(shè)為0，從而形成稀疏矩陣，通過(guò)對(duì)其進(jìn)行計(jì)算和存儲(chǔ)，達(dá)到減少內(nèi)存和計(jì)算量的目的，但是，這仍然需要先將靈敏度矩陣計(jì)算出來(lái).姚長(zhǎng)利等(2003)提出了幾何格架的概念，他們將地下劃分為規(guī)則的幾何形體，計(jì)算模型的部分幾何格架的值，隨后利用其平移等效性特征推算其他幾何格架的值.這樣，僅需要計(jì)算和存儲(chǔ)部分幾何格架的值，就可以推算出所有幾何格架的值，減少了計(jì)算量.與其類似，本文從重力和重力梯度正演公式的角度出發(fā)，提出一種快速計(jì)算靈敏度矩陣的方法.文中將正演公式看作是包含三個(gè)變量的函數(shù)式，統(tǒng)計(jì)這三個(gè)變量的值，將其依次代入公式中，求解得到構(gòu)成靈敏度矩陣的值.按照一定的規(guī)律將數(shù)值存到另一矩陣中，依據(jù)小長(zhǎng)方體的大小和位置信息，計(jì)算得到其對(duì)應(yīng)數(shù)值在矩陣中的位置，取出該值構(gòu)建靈敏度矩陣.這一方法對(duì)正演公式的變化量進(jìn)行統(tǒng)計(jì)，避免了重復(fù)計(jì)算，從而縮短了計(jì)算時(shí)間.更重要的是，我們可以將計(jì)算得到的矩陣和小長(zhǎng)方體的維度信息存儲(chǔ)在數(shù)據(jù)庫(kù)中，下次遇到相同的情況，可以直接進(jìn)行提取.

對(duì)于重力和重力梯度數(shù)據(jù)的反演問(wèn)題，文中選擇有限內(nèi)存BFGS擬牛頓法來(lái)進(jìn)行求解.這是一種可以用來(lái)解決大規(guī)模反演問(wèn)題的優(yōu)化方法，在三維電磁反演中已經(jīng)得到了廣泛的應(yīng)用(Newman and Boggs，2004；Haber，2005；Avdeeva and Avdeev，2006；Plessix and Mulder，2008；Avdeev and Avdeeva，2009；Zhang et al.，2012；劉云鶴和殷長(zhǎng)春，2013).這一方法僅需要計(jì)算目標(biāo)函數(shù)的梯度，其核心是通過(guò)幾組修正向量來(lái)構(gòu)建Hessian近似矩陣，從而降低了對(duì)內(nèi)存的要求.

為了解決重力和重力梯度聯(lián)合反演問(wèn)題，還需要注意深度加權(quán)函數(shù)的選擇.因?yàn)橹亓椭亓μ荻葦?shù)據(jù)的核函數(shù)隨深度增加是衰減的，直接對(duì)數(shù)據(jù)進(jìn)行反演計(jì)算，會(huì)產(chǎn)生趨膚效應(yīng).而且，因?yàn)橹亓椭亓μ荻葦?shù)據(jù)的核函數(shù)隨深度增加其衰減速率不同，所以，通過(guò)近似核函數(shù)衰減速率提出的深度加權(quán)函數(shù)不能同時(shí)應(yīng)用在重力和重力梯度數(shù)據(jù)上.因此，本文選擇了Commer(2011)提出的基于異常體的分布的加權(quán)函數(shù)，這一加權(quán)函數(shù)需要異常體的大致埋深信息.對(duì)于異常體之上到地表的部分，給予較小的權(quán)值，對(duì)于異常體之下的深部則給予較大的權(quán)值.為了獲得異常體的埋深信息，我們首先選擇單一分量進(jìn)行反演，隨后依據(jù)反演結(jié)果對(duì)異常體的埋深進(jìn)行評(píng)估.在本文中，我們通過(guò)再加權(quán)的方法將加權(quán)函數(shù)結(jié)合到擬合函數(shù)上.

本文首先對(duì)正演原理進(jìn)行介紹，給出快速計(jì)算靈敏度矩陣的方法，之后引入反演方法，內(nèi)容包括：穩(wěn)定函數(shù)的選擇及相應(yīng)公式的推導(dǎo)，有限內(nèi)存BFGS擬牛頓法的流程，深度加權(quán)引入，再加權(quán)優(yōu)化求解，正則化參數(shù)的選取.文中通過(guò)模型試驗(yàn)對(duì)有限內(nèi)存BFGS擬牛頓法的特點(diǎn)及正則化參數(shù)選取對(duì)反演過(guò)程的影響做了討論，并探究重力和重力梯度數(shù)據(jù)聯(lián)合反演的優(yōu)勢(shì).最后，利用美國(guó)路易斯安娜州文頓巖丘的實(shí)測(cè)重力和重力梯度數(shù)據(jù)對(duì)方法進(jìn)行驗(yàn)證.

2正演介紹

重力正演的目的是計(jì)算地下密度分布在地上引起的響應(yīng)，本文在笛卡兒坐標(biāo)系下進(jìn)行計(jì)算.假定地下存在一個(gè)長(zhǎng)方體，其長(zhǎng)寬高分別沿X、Y、Z軸方向，其中Z軸豎直向下.長(zhǎng)方體在三個(gè)方向的坐標(biāo)變化范圍是(x1,x2)，(y1,y2)，(z1,z2).觀測(cè)點(diǎn)坐標(biāo)為r=(x,y,z)，則依據(jù)Zhdanov等(2004)，可以得到重力數(shù)據(jù)的正演表達(dá)式為

(1)

其中γ是萬(wàn)有引力常數(shù)，一般取γ=6.67×10-11m3/(kg·s2)，ρ表示地下異常體的密度分布.重力梯度數(shù)據(jù)表達(dá)式為

(2)

其中

α，β=x,y,z.

(3)

為了進(jìn)一步將重力和重力梯度數(shù)據(jù)的表達(dá)式簡(jiǎn)化，文中采用了Forsberg(1984)的方法，引入了一個(gè)三元組(an,bn,cn)l，其中n=1,2，將x-xn，y-yn和z-zn循環(huán)代入可以得到式(4)—(6)：

l=1∶(an,bn,cn)1=(y-yn,z-zn,x-xn)，

(4)

l=2∶(an,bn,cn)2=(z-zn,x-xn，y-yn)，

(5)

l=3∶(an,bn,cn)3=(x-xn，y-yn，z-zn)，

(6)

其中，l=1,2,3分別對(duì)應(yīng)x,y,z.因而有g(shù)z=g3，由式(1)和式(6)得到重力數(shù)據(jù)的正演表達(dá)式為

-ailn(rijk+bj)-bjln(rijk+ai)]，

(7)

(8)

類似的，對(duì)于l,l′=1,2,3且l≠l′，可以得到其他重力梯度分量的表達(dá)式如下：

(9)

假定測(cè)區(qū)的大小為L(zhǎng)x×Ly×Lz，地下被劃分為N=nx×ny×nz個(gè)小長(zhǎng)方體，每個(gè)小長(zhǎng)方體的大小為dx×dy×dz.數(shù)據(jù)觀測(cè)點(diǎn)規(guī)則分布，其在地面的投影位于小長(zhǎng)方體的中心，觀測(cè)點(diǎn)的坐標(biāo)為(x,y,z)，觀測(cè)點(diǎn)個(gè)數(shù)為Nobs=nx×ny.假定密度ρ是單位密度，則由式(7)、(8)、(9)計(jì)算每個(gè)小長(zhǎng)方體對(duì)觀測(cè)點(diǎn)的影響，得到一個(gè)大小為Nobs×N的矩陣，用A表示，稱之為靈敏度矩陣.用d表示觀測(cè)數(shù)據(jù)，m表示地下密度分布，則正演的矩陣形式為

d=Am,

(10)

式中，d的大小為Nobs×1，m的大小為N×1.

在正演計(jì)算中，靈敏度矩陣的計(jì)算消耗了大量時(shí)間，為此文中給出了快速計(jì)算靈敏度矩陣的方法.

因?yàn)榈叵乱?guī)則劃分，所以對(duì)于(an,bn,cn)l，n=1,2，有：

(11)

式中d1,d2,d3是對(duì)應(yīng)方向上的網(wǎng)格的大小，是常數(shù).由此，式(7)、(8)、(9)包含三個(gè)變量a1,b1,c1，其統(tǒng)一的表達(dá)式為：

(12)

令a1=x-x1，b1=y-y1，c1=c-c1，x1、y1、z1是小長(zhǎng)方體中的點(diǎn)的坐標(biāo)的最小值.快速計(jì)算靈敏度矩陣的方法的流程如下：

(2) 將X、Y、Z中不同的數(shù)據(jù)進(jìn)行組合，組合種類為Nk=(2nx-1)×(2ny-1)×nz，則每個(gè)分量的數(shù)據(jù)量為Nk，將不同組合代入式(7)—(9)，得到組成重力和重力梯度數(shù)據(jù)的靈敏度矩陣的值.

(3)依據(jù)X、Y、Z中的數(shù)據(jù)計(jì)算出一個(gè)索引值，記錄數(shù)據(jù)的存放位置.對(duì)于任意分量，創(chuàng)建一個(gè)矩陣Vα，其大小為Nk×1.對(duì)于第k次計(jì)算，X、Y、Z中數(shù)據(jù)對(duì)應(yīng)的索引分別為ix，iy，iz，則數(shù)據(jù)在矩陣Vα中的位置ik由式(13)給出.

ik=(ix-1)×(2ny-1)×nz+(iy-1)×nz+iz.

(13)

(4) 依據(jù)計(jì)算得到的值來(lái)構(gòu)建靈敏度矩陣.對(duì)于任意分量，其靈敏度矩陣記為A，Aij表示標(biāo)記為j的小長(zhǎng)方體對(duì)觀測(cè)點(diǎn)i的影響.觀測(cè)點(diǎn)坐標(biāo)為(xi,yi,zi)，小長(zhǎng)方體上的點(diǎn)在三個(gè)方向坐標(biāo)取值的最小值是(xj,yj,zj)，則ix，iy，iz可用式(14)、(15)、(16)得到

(14)

(15)

iz=1+((zi-zj)-(z-(nz-1)dz))/dz.

(16)

將ix、iy和iz代入式(13)可以得到Aij的值在矩陣Vα中的位置，進(jìn)而得出Aij的值.

因?yàn)閂α中的值是按照一定規(guī)律存儲(chǔ)的，其數(shù)值與網(wǎng)格大小有關(guān)，所以可以將常用到的幾種網(wǎng)格大小得到的Vα單獨(dú)存儲(chǔ)，建立一個(gè)數(shù)據(jù)庫(kù)，下次使用時(shí)，可以直接調(diào)用.

3反演方法

3.1穩(wěn)定函數(shù)的選擇

通常，對(duì)于重力和重力梯度數(shù)據(jù)的線性反演問(wèn)題，用下面的式子表示：

(17)

(18)

這里α是正則化參數(shù)，φd是擬合函數(shù)，φs是穩(wěn)定函數(shù)，l和u代表模型參數(shù)m的取值范圍.對(duì)于正則化反演問(wèn)題，我們可以通過(guò)數(shù)值優(yōu)化算法來(lái)預(yù)測(cè)模型參數(shù)m的值，進(jìn)而達(dá)到最小化目標(biāo)函數(shù)的目的.

Last等(1983)提出了最小支撐(MS)函數(shù)，通過(guò)尋找具有最小體積的異常體來(lái)解釋重力異常.基于這一函數(shù)，Portniaguine等(1999)提出了最小梯度支撐(MGS)算子并將其應(yīng)用到了重力和可控源電磁數(shù)據(jù)反演中.MGS算子能夠用來(lái)尋找具有大的梯度的結(jié)構(gòu)體，并尋找其最小體積，對(duì)于具有陡峭邊界的地質(zhì)體，能夠獲得一個(gè)清晰、準(zhǔn)確的邊界.基于Zhdanov(2009)對(duì)該算子的改進(jìn)，我們得到穩(wěn)定函數(shù)的表達(dá)式為

(19)

φ=φd+αφs

l≤m≤u.

(20)

基于式(20)，可以給出重力和重力梯度數(shù)據(jù)五個(gè)獨(dú)立分量聯(lián)合反演的目標(biāo)函數(shù)的表達(dá)式為

(21)

3.2有限內(nèi)存擬牛頓法

有限內(nèi)存BFGS擬牛頓法(LBFGS)是一種適用于解決大規(guī)模優(yōu)化問(wèn)題的反演方法(Nocedal and Wright，1999).這一方法僅要求計(jì)算目標(biāo)函數(shù)的梯度.對(duì)于一般的擬牛頓法，比如BFGS擬牛頓法，需要存儲(chǔ)大小為N×N的Hessian矩陣H，對(duì)于有限內(nèi)存BFGS擬牛頓法，可以通過(guò)n(n取3～20之間的數(shù))組修正向量來(lái)構(gòu)建Hessian近似矩陣，需要的內(nèi)存空間為2×n×N，大大減少了對(duì)內(nèi)存的要求.這一方法已經(jīng)被用來(lái)解決三維電磁(EM)反演問(wèn)題(Newman and Boggs，2004； Haber，2005；Plessix and Mulder，2008)和一維和三維的大地電磁(MT)反演(Avdeeva and Avdeev，2006；Avdeev and Avdeeva，2009；Zhang et al.，2012).對(duì)于重力和重力梯度數(shù)據(jù)聯(lián)合反演問(wèn)題，與傳統(tǒng)的重力和重力梯度數(shù)據(jù)反演相比，數(shù)據(jù)量增大，因此本文選擇有限內(nèi)存BFGS法來(lái)求解該問(wèn)題.

si=mi+1-mi，

(22)

(23)

(24)

結(jié)合上面的信息，依據(jù)Nocedal等(1999)，文中給出了有限內(nèi)存BFGS擬牛頓法的流程：

(1) 設(shè)定模型參數(shù)m的初值m0=0及修正向量的組數(shù)n；

(5) 通過(guò)Wolfe準(zhǔn)則，搜索合適的步長(zhǎng)βk，計(jì)算得到mk+1=mk+βkpk；

(6) 對(duì)目標(biāo)函數(shù)的值和迭代次數(shù)進(jìn)行判斷，如果目標(biāo)函數(shù)的值收斂或者達(dá)到迭代次數(shù)則終止循環(huán)，否則轉(zhuǎn)到(2).

在有限內(nèi)存LBFGS擬牛頓法中，為了保證其收斂，文中選擇精確線搜索方法Wolfe準(zhǔn)則，其表達(dá)式為

(25)

其中0

同時(shí)，為了保證Hessian近似矩陣Hk是正定矩陣，向量sk和yk必須滿足下面的條件

(26)

(27)

3.3深度加權(quán)

為了抵消趨膚效應(yīng)，需要采用深度加權(quán)函數(shù).依據(jù)Li等(1996)，文中給出模型試驗(yàn)說(shuō)明本文的加權(quán)方式.圖1所示為一個(gè)長(zhǎng)方體模型，大小為500 m×500 m×500 m，頂面埋深為500 m，剩余密度為1 g·cm-3，為了與實(shí)際數(shù)據(jù)吻合，加入了10%的高斯噪聲.

對(duì)于三維重力和重力梯度正演模型，其表達(dá)式的矩陣形式為式(10).由式(1)和(2)，重力和重力梯度數(shù)據(jù)的核函數(shù)隨深度增加是衰減的，因此會(huì)引起反演結(jié)果的趨膚效應(yīng).文中引入加權(quán)函數(shù)矩陣W來(lái)抵消這一影響.式(10)轉(zhuǎn)換為

(28)

同樣式(17)轉(zhuǎn)換為

(29)

(30)

對(duì)于重力數(shù)據(jù)，取β=1，對(duì)w(z)進(jìn)行歸一化處理，得到w(z)隨深度變化如圖2中的實(shí)線所示，將其應(yīng)用到gz分量的反演計(jì)算中，得到結(jié)果如圖1b所示，圖中顯示反演結(jié)果克服了趨膚效應(yīng).對(duì)于重力梯度數(shù)據(jù)，取β=1.5，對(duì)w(z)進(jìn)行歸一化處理，得到w(z)隨深度變化如圖2中虛點(diǎn)線所示.將其應(yīng)用到gzz分量的反演計(jì)算中，得到結(jié)果如圖1c所示，圖中顯示趨膚效應(yīng)被克服了.

圖1　單個(gè)長(zhǎng)方體模型和不同分量的反演結(jié)果的剖面圖(a) 模型三維立體圖； (b) gz分量反演結(jié)果； (c) gzz分量反演結(jié)果； (d) z1=400時(shí)， gz和gzz聯(lián)合反演結(jié)果； (e) z1=300時(shí)，gz和gzz聯(lián)合反演結(jié)果； (f) z1=500時(shí)， gz和gzz聯(lián)合反演結(jié)果.Fig.1　Single prism model and profiles of inversion results of different components(a) 3D model； (b) Inversion result of gz component； (c) Inversion result of gzz component； (d) z1=400， inversion result of gzand gzz component； (e) z1=300，inversion result of gz and gzz component； (f) z1=500，inversion result of gz and gzz component.

圖2　不同深度加權(quán)函數(shù)隨深度變化曲線Fig.2　Weighting functions varying with depth

由圖1b和圖1c得出，對(duì)目標(biāo)函數(shù)的擬合函數(shù)部分進(jìn)行加權(quán)消除趨膚效應(yīng)是可行的.因?yàn)橹亓?shù)據(jù)和重力梯度數(shù)據(jù)具有不同的衰減速率，所以，從理論上來(lái)說(shuō)，通過(guò)近似衰減速率得到的加權(quán)函數(shù)不能直接應(yīng)用到重力和重力梯度數(shù)據(jù)聯(lián)合反演中.Commer等(2011)提出了一種空間加權(quán)函數(shù)，Commer(2011)將其應(yīng)用到了重力數(shù)據(jù)的反演計(jì)算中，這一函數(shù)需要了解模型的大致埋深，隨后針對(duì)異常在垂直方向的分布規(guī)律，對(duì)于近地表的部分，給予較小的權(quán)值以達(dá)到阻尼的效果，對(duì)于可能存在異常的區(qū)域，則給予較大的權(quán)值.因?yàn)樵诶硐霔l件下，對(duì)同一個(gè)地質(zhì)體的重力和重力梯度數(shù)據(jù)分別進(jìn)行反演計(jì)算，得到的模型是一致的.所以，這一基于模型分布提出的加權(quán)函數(shù)，可以用于重力和重力梯度數(shù)據(jù)的聯(lián)合反演.其表達(dá)式為：

(31)

其中α和r是常數(shù)，z1是與深度信息有關(guān)的常數(shù)，文中取α=0.001，r=1.其中α的值由Commer(2011)給出，是經(jīng)驗(yàn)值.依據(jù)劉銀萍等(2013)，α也可以取

其他很小的值，它的大小直接決定了近地表頂層壓制作用的大小，其值越小，壓制作用越大，反演結(jié)果越集中，反之亦然.對(duì)于參數(shù)r，其主要作用是保證當(dāng)z→0時(shí)，w(z)的值是一個(gè)很小的值或者保證z→0時(shí)，w(z)≈α，從而影響近地表頂層壓制作用的大小.由圖1b和圖1c，取z1=400，得到深度加權(quán)函數(shù)如圖2中虛線所示，將其應(yīng)用到gz和gzz分量的聯(lián)合反演中，得到圖1d.可以看到，反演結(jié)果克服了趨膚效應(yīng)，能夠準(zhǔn)確反映出異常體的位置.

由此發(fā)現(xiàn)，通過(guò)式(29)對(duì)單一分量進(jìn)行反演計(jì)算得到的異常體的深度值小于實(shí)際深度.這是因?yàn)闆](méi)有其他約束條件，得到的反演結(jié)果(圖1b和圖1c)比較發(fā)散.但是將這一信息結(jié)合到式(31)中，將其作為深度加權(quán)函數(shù)應(yīng)用到反演計(jì)算，得到的反演結(jié)果與真實(shí)模型一致.

為了進(jìn)一步驗(yàn)證預(yù)測(cè)深度值z(mì)1對(duì)反演結(jié)果的影響，分別取z1=300和z1=500，得到反演結(jié)果如圖1e和圖1f.由圖1e，當(dāng)z1=300時(shí)，異常體的位置略有上升，但是其高密度部分(>0.25 g·cm-3)與真實(shí)模型一致.當(dāng)z1=500時(shí)，深度信息是準(zhǔn)確的，圖1f中異常體的位置與真實(shí)模型一致.由此，雖然單一分量反演得到的異常體的深度信息與異常體實(shí)際埋深存在一定的差異，但是，將這一信息作為先驗(yàn)信息結(jié)合到深度加權(quán)函數(shù)中，進(jìn)行實(shí)際的反演計(jì)算，得到的反演結(jié)果是可信的.即使預(yù)測(cè)的深度在一定范圍內(nèi)波動(dòng)，也能得到可信的結(jié)果.同時(shí)，對(duì)比圖1d、圖1e和圖1f，可以發(fā)現(xiàn)，采用預(yù)測(cè)得到的深度值作為加權(quán)函數(shù)，得到的反演結(jié)果的密度的最大值更接近真實(shí)模型.

3.4再加權(quán)優(yōu)化求解

將加權(quán)函數(shù)應(yīng)用到目標(biāo)函數(shù)中，則式(20)變?yōu)?/p>

(32)

對(duì)于擬合函數(shù)部分，對(duì)mw求梯度得到

(33)

(34)

其中

(35)

(36)

對(duì)等式

(37)

進(jìn)行積分計(jì)算，并重新排列得到

+?? VCδmw·nds，

(38)

其中n是單位向量，指向區(qū)域V的外部，?V是V的表面.依據(jù)齊次諾依曼邊界條件

(39)

可以得到

(40)

則由式(34)、式(38)和式(40)得到

(41)

由式(33)和式(41)，得到目標(biāo)函數(shù)的梯度的表達(dá)式為

(42)

由式(20)，模型參數(shù)m存在邊界值，為了增強(qiáng)邊界約束，將模型參數(shù)的邊界信息加入到反演過(guò)程中，文中采用懲罰函數(shù)的方法，在每次迭代結(jié)束后，由mw計(jì)算模型參數(shù)m，對(duì)于超過(guò)邊界值的部分，將邊界值賦給相應(yīng)的位置.

3.5正則化參數(shù)的選取

在反演計(jì)算過(guò)程中，需要確定正則化參數(shù)的值.Farquharson等(2004)對(duì)依據(jù)偏差原理、GCV準(zhǔn)則和L曲線準(zhǔn)則預(yù)測(cè)正則化參數(shù)的方法進(jìn)行了對(duì)比.他們對(duì)同一數(shù)據(jù)加不同的噪聲，在保證其他反演參數(shù)不變的情況下，將這三種方法應(yīng)用到具有不同噪聲水平的數(shù)據(jù)，結(jié)果顯示依據(jù)GCV準(zhǔn)則和L曲線準(zhǔn)則獲得的結(jié)果比依據(jù)偏差原理獲得的結(jié)果更準(zhǔn)確.進(jìn)而得出，依據(jù)GCV準(zhǔn)則和L曲線準(zhǔn)則預(yù)測(cè)正則化參數(shù)的方法能夠在反演計(jì)算中給出恰當(dāng)?shù)恼齽t化參數(shù)的值.Portniaguine等(2002)給出了一種簡(jiǎn)單的數(shù)值方法，給定正則化參數(shù)的初值，在迭代計(jì)算中，正則化參數(shù)的值依次遞減.Farquharson等(2004)給出的反演計(jì)算過(guò)程中的正則化參數(shù)的值，也符合這一規(guī)律.依據(jù)Zhdanov(2002)，目標(biāo)函數(shù)是關(guān)于正則化參數(shù)的單調(diào)函數(shù)，正則化參數(shù)的值減小，目標(biāo)函數(shù)的值也減小.我們選擇數(shù)值優(yōu)化算法來(lái)求解目標(biāo)函數(shù)時(shí)，一般是通過(guò)最小化目標(biāo)函數(shù)來(lái)尋找最優(yōu)解.對(duì)于每次迭代獲得的解，如果正則化參數(shù)的值減小，則相應(yīng)的目標(biāo)函數(shù)的值也減小，加快了收斂速度.

下面，在模型試驗(yàn)中，在保持其他參數(shù)設(shè)置不變的情況下，分別將α取1和Portniaguine等(2002)提出的正則化參數(shù)的計(jì)算方法應(yīng)用到算法中，對(duì)含不同噪聲水平的數(shù)據(jù)進(jìn)行反演計(jì)算.Portniaguine等(2002)提出的正則化參數(shù)的計(jì)算方法如式(43)和式(44)所示.

(43)

(44)

其中，擬合函數(shù)φd和穩(wěn)定函數(shù)φs的表達(dá)式可以由式(32)得到，初次迭代計(jì)算后得到的正則化參數(shù)的表達(dá)式為式(43)，隨后的迭代計(jì)算中，則通過(guò)式(44)計(jì)算正則化參數(shù)的值.

4模型試驗(yàn)

首先通過(guò)單一模型試驗(yàn)來(lái)說(shuō)明有限內(nèi)存BFGS擬牛頓法的特點(diǎn)，隨后通過(guò)組合模型試驗(yàn)1來(lái)比較正則化參數(shù)取值方式不同對(duì)反演結(jié)果的影響，并討論不同噪聲水平數(shù)據(jù)對(duì)正則化參數(shù)選取的影響，最后通過(guò)組合模型試驗(yàn)2來(lái)探究重力和重力梯度數(shù)據(jù)聯(lián)合反演的優(yōu)勢(shì).

4.1單一模型試驗(yàn)

建立如下的模型：如圖3所示，模型由一個(gè)長(zhǎng)方體組成，其大小為4000 m×4000 m×300 m.長(zhǎng)方體的頂面埋深為400 m，剩余密度為1 g·cm-3.測(cè)區(qū)的大小為10000 m×10000 m×1000 m，地下被劃分為100×100×10個(gè)小立方體，每個(gè)小立方體的大小為100 m×100 m×100 m，觀測(cè)點(diǎn)的個(gè)數(shù)為100×100個(gè).為了更貼近實(shí)際測(cè)量，加入占最大值10%的高斯噪聲.圖4a和圖4b分別給出了觀測(cè)數(shù)據(jù)gz和gzz的等值線圖.

圖3　模型三維透視圖，其大小為4000 m×4000 m×300 m,頂面埋深為400 m,其剩余密度為1 g·cm-3Fig.3　Perspective of the model. The size of the model is 4000 m×4000 m×300 m. The depth to the top of the model is 400 m.The density contrast of the model is 1 g·cm-3

選擇gzz數(shù)據(jù)，結(jié)合式(29)進(jìn)行反演計(jì)算，其中

加權(quán)函數(shù)采用式(30)，β=1.5，得到反演結(jié)果在x=5000 m處的剖面如圖5b所示，圖5a為理論模型在x=5000 m處的剖面.由圖5b，異常分布在400 m以下的概率很大，這一信息與理論模型相符.由此得到式(31)的參數(shù)設(shè)置為α=0.001，r=1，z1=400.

圖6a—6c分別給出了gz|gzz、gz|gxz|gyz|gzz和gz|gxy|gxz|gyy|gyz|gzz三組分量組合的反演結(jié)果，可以看出，三組分量的反演結(jié)果類似，與圖5a對(duì)比，圖6a—6c結(jié)果與理論模型相吻合.

此處的模型試驗(yàn)中，單一分量的測(cè)點(diǎn)個(gè)數(shù)為10000，網(wǎng)格數(shù)為100000.在gz|gzz反演中，計(jì)算用到的測(cè)點(diǎn)數(shù)據(jù)個(gè)數(shù)為20000.相應(yīng)的，在gz|gxz|gyz|gzz和gz|gxy|gxz|gyy|gyz|gzz反演中，計(jì)算用到的測(cè)點(diǎn)數(shù)據(jù)個(gè)數(shù)分別為40000和60000.對(duì)于gz|gzz，反演計(jì)算中迭代過(guò)程需要468.196 s，對(duì)于gz|gxz|gyz|gzz，反演計(jì)算中迭代過(guò)程需要694.870 s，對(duì)于gz|gxy|gxz|gyy|gyz|gzz，反演計(jì)算中迭代過(guò)程需要1623.908 s.對(duì)于圖3中的模型，每一分量的靈敏度矩陣大小為10000×100000，存儲(chǔ)單一分量靈敏度矩陣需要存儲(chǔ)空間為5.36 GB，采用文中提出的方法，存儲(chǔ)組成靈敏度矩陣的值，需要的存儲(chǔ)空間為1.6 MB.對(duì)于gz|gzz反演，式(21)中矩陣A的大小為20000×100000，對(duì)于gz|gxz|gyz|gzz反演，A的大小為40000×100000，對(duì)于gz|gxy|gxz|gyy|gyz|gzz反演， A的大小為60000×100000.與單一分量反演相比較，多分量組合反演的數(shù)據(jù)量和數(shù)據(jù)規(guī)模都大大增加，因此文中選擇有限內(nèi)存BFGS擬牛頓法對(duì)反演問(wèn)題進(jìn)行求解.

圖4　(a) gz的等值線圖，包含最大值10%的高斯噪聲； (b) gzz的等值線圖，包含最大值10%的高斯噪聲Fig.4　(a) Contour map of gz, contaminated by 10% Gaussian noise； (b) Contour map of gzz, contaminated by 10% Gaussian noise

圖5　(a) 模型在x=5000 m處的垂直剖面圖， (b) gzz反演結(jié)果在x=5000 m處的垂直剖面圖，結(jié)果顯示異常體在地下400 mFig.5　(a) The vertical section of the synthetic model at x=5000 m， (b) the vertical section of gzz inversion result at x=5000 m, it shows that the anomalous body is located 400 m below the surface

圖6　(a) gz|gzz反演結(jié)果的垂直剖面，(b) gz|gxz|gyz|gzz反演結(jié)果的垂直剖面，(c) gz|gxy|gxz|gyy|gyz|gzz反演結(jié)果的垂直剖面，垂直剖面的位置為x=5000 mFig.6　(a) The vertical section of gz|gzz inversion result, (b) the vertical section of gz|gxz|gyz|gzz inversion result, (c) the vertical section of gz|gxy|gxz|gyy|gyz|gzz inversion result. The location of the vertical section is x=5000 m

由3.2節(jié)中有限內(nèi)存BFGS擬牛頓法的流程可知，在求解反演問(wèn)題的過(guò)程中，使用多組修正向量來(lái)構(gòu)建Hessian矩陣是該方法的關(guān)鍵點(diǎn).Avdeev等(2009)指出有限內(nèi)存BFGS擬牛頓法用n組修正向量來(lái)構(gòu)建Hessian矩陣，它對(duì)內(nèi)存的要求與2×n×N成正比，其中N代表網(wǎng)格數(shù)，這要求遠(yuǎn)遠(yuǎn)小于其他類型牛頓法的要求，其他類型的牛頓法對(duì)內(nèi)存要求正比于N×M或N×N，其中M代表測(cè)點(diǎn)的個(gè)數(shù).對(duì)圖3中模型，Hessian矩陣的大小為100000×100000，如果對(duì)這一矩陣進(jìn)行存儲(chǔ)，需要的空間大于30 GB.通過(guò)有限內(nèi)存BFGS擬牛頓法，我們僅需要3.64 MB的內(nèi)存空間，這里n=5.由此，與其他類型的牛頓法相比，有限內(nèi)存BFGS擬牛頓法節(jié)省了大量?jī)?nèi)存空間.

4.2組合模型試驗(yàn)1

Farquharson等(2004)通過(guò)二維模型和三維模型對(duì)反演過(guò)程中正則化參數(shù)計(jì)算進(jìn)行了論述，這里，建立類似的三維模型：如圖7a，模型由兩個(gè)長(zhǎng)方體組成，其中紅色長(zhǎng)方體大小為400 m×200 m×200 m,頂面埋深為400 m,剩余密度為1 g·cm-3，橘紅色長(zhǎng)方體大小為400 m×200 m×400 m，頂面埋深為400 m,剩余密度為0.8 g·cm-3.測(cè)區(qū)大小為1200 m×1200 m×1200 m,將地下劃分為12×12×12個(gè)小長(zhǎng)方體，每個(gè)小長(zhǎng)方體大小為100 m×100 m×100 m,測(cè)點(diǎn)個(gè)數(shù)為12×12.為了觀察不同噪聲水平數(shù)據(jù)對(duì)正則化參數(shù)的影響，將最大值3%,5%和10%的高斯噪聲分別加入到觀測(cè)數(shù)據(jù)中.圖7b,7c和7d分別給出了含3%,5%和10%高斯噪聲的gzz的等值線圖.

因?yàn)樵?.1節(jié)中，gz|gzz分量組合能夠得到清晰的反演結(jié)果且計(jì)算時(shí)間短，所以這里選擇gz|gzz分量組合進(jìn)行反演計(jì)算.加權(quán)函數(shù)參數(shù)取α=0.001，r=1，z1=400.

圖8a，8b分別給出了理論模型在z=500 m的橫切面和x=500 m處的縱切面.圖8c—8h中的反演結(jié)果是在正則化參數(shù)取1的情況下，分別對(duì)含3%、5%和10%噪聲水平的數(shù)據(jù)進(jìn)行反演計(jì)算得到的.對(duì)比發(fā)現(xiàn)，三組數(shù)據(jù)反演得到的結(jié)果相近，與真實(shí)模型吻合.由此說(shuō)明對(duì)不同噪聲水平數(shù)據(jù)，正則化參數(shù)取1可行.通過(guò)式(43)和式(44)計(jì)算正則化參數(shù)，隨后分別對(duì)三組數(shù)據(jù)進(jìn)行反演計(jì)算，其他參數(shù)設(shè)置保持不變，得到的反演結(jié)果如圖9a—9f所示.與圖8c—8h比較，圖9中反演結(jié)果與圖8中反演結(jié)果一致，且圖9中的反演結(jié)果與實(shí)際模型吻合程度更好.由此說(shuō)明，在實(shí)際計(jì)算中，依據(jù)式(43)和(44)計(jì)算正則化參數(shù)，能夠得到更好的反演結(jié)果.

文中對(duì)正則化參數(shù)取值方式不同時(shí)，反演算法的收斂性進(jìn)行比較.由圖10中的結(jié)果可以看出，當(dāng)正則化參數(shù)隨迭代過(guò)程變化時(shí)，目標(biāo)函數(shù)的收斂速度更快，數(shù)據(jù)能夠更好的擬合.

圖7　(a)模型的三維透視圖，其中紅色模型密度為1 g·cm-3，橘紅色模型密度為0.8 g·cm-3； (b) gzz的等值線圖，包含最大值3%的噪聲； (c) gzz的等值線圖，包含最大值5%的噪聲； (d) gzz的等值線圖，包含最大值10%的噪聲Fig.7　(a) Perspective of the model. The red blocks have a density of 1 g·cm-3, the darker blocks have a density of 0.8 g·cm-3， (b) Contour map of gzz, contaminated by 3% Gaussian noise; (c) Contour map of gzz, contaminated by 5% Gaussian noise; (d) Contour map of gzz, contaminated by 10% Gaussian noise

圖8　理論模型的橫切面(a)和垂直切面(b),正則化參數(shù)取1， gz|gzz反演結(jié)果的橫切面(c，e，g)，gz|gzz反演結(jié)果的垂直切面(d，f，h).左邊一列橫切面的深度為z=500 m，右邊一列垂直切面的位置為x=500 mFig.8　Depth section of the true model (a), vertical section of the true model (b), regularization parameter is set equal to 1, the depth section of the gz|gzz inversion (c,e,g), and vertical section of the gz|gzz inversion (d,f,h). The left column is at z=500 m. The right column is at x=500 m.

圖9　依據(jù)式(43)和式(44)計(jì)算正則化參數(shù)，gz|gzz反演結(jié)果的橫切面(a，c，e)，gz|gzz反演結(jié)果的垂直切面(b，d，f).左邊一列橫切面的深度為z=500 m，右邊一列垂直切面的位置為x=500 mFig.9　Regularization parameter was calculated based on equation (43) and (44). The depth section of gz|gzz incersion (a,c,e), the vertical section of gz|gzz inversion (b,d,f). The left column is at z=500 m, the right column is at x=500 m

圖10　對(duì)包含不同噪聲(a)3%， (b) 5%， (c) 10%的gz|gzz組合進(jìn)行反演計(jì)算，得到的反演過(guò)程中目標(biāo)函數(shù)的變化實(shí)線表示正則化參數(shù)取1時(shí)的計(jì)算結(jié)果，帶點(diǎn)的實(shí)線表示正則化參數(shù)變化時(shí)的計(jì)算結(jié)果.Fig.10　Values of the objective function during inversion of gz|gzz components. The noise level is (a) 3%, (b) 5%, (c) 10% of the maximum valueHere, we set the regularization parameter equal to 1, the result is drawn with a solid line. Calculated the regularization parameter with equation (43) and (44), the result is drawn in solid line with points.

圖11　(a)當(dāng)正則化參數(shù)取1時(shí)，不同噪聲水平的gz|gzz數(shù)據(jù)在反演過(guò)程中，目標(biāo)函數(shù)的變化；(b)依據(jù)式(43)和(44)計(jì)算正則化參數(shù)時(shí)，不同噪聲水平的gz|gzz數(shù)據(jù)在反演過(guò)程中，目標(biāo)函數(shù)的變化； (c)當(dāng)正則化參數(shù)變化時(shí)，在反演過(guò)程中計(jì)算得到的正則化參數(shù)圖中實(shí)線表示3%噪聲，帶點(diǎn)的實(shí)線表示5%噪聲，帶“+”的實(shí)線表示10%噪聲.Fig.11　(a) Setting the regularization parameter equal to 1, the value of the objective function for gz|gzz inversion with different noise. (b) Based on equation (43) and (44) to calculate the regularization parameter, the objective function during the gz|gzz inversion. (c) The regularization parameter during the inversionHere, the solid line represent 3% noise, the solid line with points represent 5% noise, the solid line with plus sign represent 10% noise.

下面，對(duì)比不同噪聲水平數(shù)據(jù)對(duì)反演算法的收斂性影響.由圖11a，11b可知，對(duì)不同噪聲水平數(shù)據(jù)進(jìn)行反演計(jì)算，目標(biāo)函數(shù)的收斂速度不同.圖11c則給出了正則化參數(shù)變化的情況下，對(duì)三組數(shù)據(jù)分別進(jìn)行反演計(jì)算，得到的迭代過(guò)程中的正則化參數(shù).可以看出，對(duì)于不同噪聲水平的數(shù)據(jù)，正則化參數(shù)的值在每次迭代中是相近的.由此發(fā)現(xiàn)，不同噪聲水平數(shù)據(jù)對(duì)正則化參數(shù)的取值影響很小.

基于上面的討論，可以發(fā)現(xiàn)，雖然依據(jù)式(43)和(44)，在迭代過(guò)程中更新正則化參數(shù)，能夠得到更好的反演結(jié)果，但是正則化參數(shù)取1時(shí)也能得出合理的反演結(jié)果.同時(shí)，基于3.5節(jié)中的分析，為了保持方程在迭代過(guò)程中的一致性，文中取正則化參數(shù)值為1.4.3組合模型試驗(yàn)2

建立如下模型：如圖12所示，模型由五個(gè)長(zhǎng)方體組成，其中兩個(gè)紅色長(zhǎng)方體大小為400 m×400 m×400 m，另一個(gè)紅色長(zhǎng)方體大小為200 m×200 m×200 m，它們的頂面埋深為400 m，剩余密度1 g·cm-3，兩個(gè)橘紅色長(zhǎng)方體大小分別為1000 m×100 m×500 m和800 m×200 m×200 m，它們的頂面埋深為400 m，剩余密度為0.8 g·cm-3.測(cè)區(qū)的大小為2500 m×2500 m×1500 m.將地下劃分為25×25×15個(gè)小長(zhǎng)方體，每個(gè)小長(zhǎng)方體大小為100×100×100，觀測(cè)點(diǎn)的個(gè)數(shù)為25×25，為了更貼近實(shí)際測(cè)量數(shù)據(jù)，加入了占最大值10%的高斯噪聲.圖13a和圖13b分別給出了觀測(cè)數(shù)據(jù)gz和gzz的等值線圖.

選擇gzz數(shù)據(jù)，結(jié)合式(29)進(jìn)行反演計(jì)算，其中加權(quán)函數(shù)采用式(30)，β=1.5，得到反演結(jié)果在x=900 m和x=1800 m處的垂直剖面圖，如圖14c和圖14d所示.圖14a和圖14b給出了模型在相同位置的垂直剖面圖.由圖14c，異常體分布在400 m以下的概率很大，這一信息與實(shí)際模型分布一致，圖14d中的結(jié)果顯示，異常分布在900 m以下.由圖13a和圖13b中異常值的分布，選擇圖14c中的結(jié)果作為預(yù)測(cè)的異常體的深度信息，則加權(quán)函數(shù)參數(shù)取α=0.001，r=1，z1=400.將加權(quán)函數(shù)應(yīng)用到不同重力和重力梯度分量組合的聯(lián)合反演中，得到結(jié)果如圖15和圖16所示.

圖15a—15c分別給出了理論模型在z=500 m，z=600 m和z=800 m處的橫切面圖.圖15d—15f給出了gz反演結(jié)果，從圖中，無(wú)法識(shí)別出不同的異常體.圖15g—15i給出了gzz反演結(jié)果，從圖15g中，可以識(shí)別出五個(gè)不同的異常體，其位置與真實(shí)模型的位置一致.對(duì)gz和gzz進(jìn)行聯(lián)合反演，結(jié)果如圖15j—15l所示，gz|gzz反演的結(jié)果與gzz反演結(jié)果類似.表1中給出了高斯噪聲的標(biāo)準(zhǔn)差和由不同反演結(jié)果得到的各分量的預(yù)測(cè)數(shù)據(jù)和觀測(cè)數(shù)據(jù)的差的標(biāo)準(zhǔn)差.可以發(fā)現(xiàn)gz|gzz聯(lián)合反演得到的各分量預(yù)測(cè)數(shù)據(jù)和觀測(cè)數(shù)據(jù)的差的標(biāo)準(zhǔn)差(gxx：2.737，gxy：0.887，gxz：2.440，gyy：1.907，gyz：2.259，gzz：3.671，gz：0.177)大于gzz反演結(jié)果得到的各分量預(yù)測(cè)數(shù)據(jù)和觀測(cè)數(shù)據(jù)的差的標(biāo)準(zhǔn)差(gxx：2.731，gxy：0.883，gxz：2.438，gyy：1.899，gyz：2.257，gzz：3.651，gz：0.177).這是因?yàn)間z分量的加入，識(shí)別出了gzz中與gz中包含信息不符的噪聲成分.雖然，gz|gzz聯(lián)合反演對(duì)模型的改進(jìn)很小，但是gz|gzz聯(lián)合反演得到的模型更加合理.

圖16a—16c給出了gxy|gyz|gzz反演結(jié)果的橫切面，圖16d—16f給出了gz|gxz|gyz|gzz反演結(jié)果的橫切面，這兩組反演結(jié)果類似，且與理論模型相符.將圖16d與圖15j對(duì)比，可以發(fā)現(xiàn)，圖16d中的模型更加緊致，說(shuō)明隨著反演過(guò)程中數(shù)據(jù)分量的增加，反演模型得到了改進(jìn).

圖16g—16i和圖16j—16l分別給出了gxy|gxz|gyy|gyz|gzz和gz|gxy|gxz|gyy|gyz|gzz反演結(jié)果.這兩組分量組合的反演結(jié)果與圖16a—16f類似.但是，由表1可以得出，gz|gxy|gxz|gyy|gyz|gzz反演結(jié)果得到的預(yù)測(cè)數(shù)據(jù)和觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)差和高斯噪聲的標(biāo)準(zhǔn)差最接近，說(shuō)明gz|gxy|gxz|gyy|gyz|gzz反演得到的結(jié)果更加合理.

圖12　模型三維透視圖，其中紅色模型密度為1 g·cm-3,橘紅色模型密度值為0.8 g·cm-3Fig.12　Perspective of the model. The red blocks have a density of 1 g·cm-3,the darker blocks have a density of 0.8 g·cm-3

圖13　(a) gz的等值線圖，包含最大值10%的高斯噪聲； (b) gzz的等值線圖，包含最大值10%的高斯噪聲Fig.13　(a) Contour map of gz, contaminated by 10% Gaussian noise； (b) Contour map of gzz, contaminated by 10% Gaussian noise

圖14　(a)模型在x=900 m處的垂直剖面圖； (b) 模型在x=1800 m處的垂直剖面圖； (c) gzz反演結(jié)果在x=900 m處的垂直剖面圖，結(jié)果顯示異常體在地下400 m; (d) gzz反演結(jié)果在x=1800 m處的垂直剖面圖Fig.14　(a) Profile of the synthetic model at x=900 m; (b) Profile of the synthetic model at x=1800 m; (c) Profile of gzz inversion result at x= 900 m, it shows that the anomalous body is located 450m below the surface; (d) Profile of gzz inversion result at x= 1800 m

表1　不同分量的高斯噪聲的標(biāo)準(zhǔn)差和由各分量組合得到的反演結(jié)果的預(yù)測(cè)數(shù)據(jù)和觀測(cè)數(shù)據(jù)的差的標(biāo)準(zhǔn)差　(gz: mGal，gxx|gxy|gxz|gyy|gyz|gzz:E)

基于模型試驗(yàn)，可以得出，將重力和重力梯度數(shù)據(jù)聯(lián)合反演，反演結(jié)果有一定程度的改進(jìn)，雖然改進(jìn)不大，但聯(lián)合反演得到的結(jié)果更加合理.

圖17a—17g給出了不同分量組合反演過(guò)程中目標(biāo)函數(shù)的變化，由圖中可以得出，采用有限內(nèi)存BFGS擬牛頓法最小化目標(biāo)函數(shù)，能夠?qū)崿F(xiàn)較快的收斂.

5實(shí)際數(shù)據(jù)

為了進(jìn)一步說(shuō)明方法的有效性，并說(shuō)明重力和重力梯度數(shù)據(jù)聯(lián)合反演在實(shí)際中的應(yīng)用效果，文中將方法應(yīng)用到美國(guó)路易斯安那州Vinton鹽丘的重力數(shù)據(jù)和全張量重力梯度數(shù)據(jù).Vinton鹽丘是一個(gè)刺穿型鹽丘，鹽丘上覆一個(gè)巨大的巖蓋，主要由石膏和硬石膏組成(Coker et al.，2007)，其密度約為2.75 g·cm-3(Ennen，2012).鹽丘周圍主要是沉積巖層，密度約為2.2 g·cm-3.

圖15　理論模型的橫切面(a，b，c)，gz反演結(jié)果的橫切面(d，e，f)，gzz反演結(jié)果的橫切面(g，h，i)，gz|gzz反演結(jié)果的橫切面(j，k，l).左邊一列橫切面深度為z=500 m，中間一列橫切面的深度為z=600 m，右邊一列橫切面的深度為z=800 mFig.15　Top row is depth sections of true model (a,b,c), depth sections of gz inversion (d,e,f), depth sections of gzz inversion (g,h,i), and depth sections of gz|gzz inversion. The left column is at z=500 m, the central column is at z=600 m, the right column is at z=800 m

圖16　gxz|gyz|gzz反演結(jié)果的橫切面(a，b，c)，gz|gxz|gyz|gzz反演結(jié)果的橫切面(d，e，f)，gxy|gxz|gyy|gyz|gzz反演結(jié)果的橫切面(g，h，i)，gz|gxy|gxz|gyy|gyz|gzz反演結(jié)果的橫切面(j，k，l)Fig.16　Depth sections of gxz|gyz|gzz inversion result (a,b,c), the depth section of gz|gxz|gyz|gzz inversion result (d,e,f), depth sections of gxy|gxz|gyy|gyz|gzz inversion results (g,h,i), depth sections of gz|gxy|gxz|gyy|gyz|gzz inversion result (j,k,l). The left column is at z=500 m, the centre column is at z=600 m, the right column is at z=800 m

圖17　不同分量組合反演過(guò)程中目標(biāo)函數(shù)的變化曲線：(a) gz， (b) gzz， (c) gz|gzz， (d) gxz|gyz|gzz， (e) gz|gxz|gyz|gzz， (f) gxy|gxz|gyy|gyz|gzz, (g)gz|gxy|gxz|gyy|gyz|gzzFig.17　The objective function during the inversion of different components： (a)gz， (b) gzz， (c) gz|gzz， (d) gxz|gyz|gzz， (e) gz|gxz|gyz|gzz， (f) gxy|gxz|gyy|gyz|gzz， (g) gz|gxy|gxz|gyy|gyz|gzz

數(shù)據(jù)由Bell Geospace公司在2008年測(cè)的，其中，重力數(shù)據(jù)在地面測(cè)得，全張量重力梯度數(shù)據(jù)則由FTG系統(tǒng)在飛機(jī)上測(cè)得.圖18a—18g展示了文中所用到的數(shù)據(jù)，其中圖18f為重力數(shù)據(jù)，圖18a—18e和圖18g是重力梯度數(shù)據(jù).重力梯度數(shù)據(jù)經(jīng)過(guò)地形改正獲得，地形改正的值為2.2 g·cm-3(Oliveira and Barbosa，2013).參考Geng等(2014)，我們對(duì)數(shù)據(jù)做了一個(gè)帶通濾波，達(dá)到增強(qiáng)數(shù)據(jù)信號(hào)的目的.對(duì)于重力數(shù)據(jù)，通過(guò)高通濾波去除背景場(chǎng).選定了大小為4.2 km×4.2 km的測(cè)區(qū)，將地下劃分為42×42×20個(gè)小長(zhǎng)方體，每個(gè)長(zhǎng)方體的大小為100×100×50.采用依次計(jì)算每個(gè)小長(zhǎng)方體對(duì)觀測(cè)點(diǎn)影響的方法計(jì)算gzz分量的靈敏度矩陣，需要8605.7 s，而通過(guò)本文提出的方法，僅需要344.3 s.

在模型試驗(yàn)中，首先對(duì)單一分量進(jìn)行反演，由反演結(jié)果可以得到模型的深度信息，隨后將這一信息結(jié)合到深度加權(quán)函數(shù)中，對(duì)實(shí)際數(shù)據(jù)，我們采用同樣的步驟.首先，采用gzz數(shù)據(jù)，利用式(29)和(30)進(jìn)行反演計(jì)算，得到結(jié)果如圖19，由圖19可以判斷出異常體分布在200 m以下.這一信息與異常體的實(shí)際埋深有一定的差異，但是由3.3節(jié)的模型試驗(yàn)，我們發(fā)現(xiàn)，依據(jù)式(29)和(30)得到的異常體埋深的估計(jì)值，雖然與實(shí)際埋深有一定的差異，但是，將這一信息作為參數(shù)結(jié)合到深度加權(quán)函數(shù)中，得到的結(jié)果與理論模型一致.而且，當(dāng)估計(jì)值在一定范圍內(nèi)波動(dòng)時(shí)，將其應(yīng)用到反演計(jì)算中，仍能獲得可信的解.因此，這一深度信息可以應(yīng)用于實(shí)際數(shù)據(jù)的計(jì)算.同時(shí)，這一信息與已知的地質(zhì)信息相符，所以此處深度加權(quán)函數(shù)式(31)的參數(shù)取值為α=0.001，r=1，z1=200.圖20和圖21給出了不同分量組合反演得到的結(jié)果.剖面線的位置如圖18g中黑色實(shí)線所示.

圖20a—20c展示了gz反演結(jié)果.表2列出了不同分量組合得到的反演結(jié)果的預(yù)測(cè)數(shù)據(jù)和觀測(cè)數(shù)據(jù)的差的預(yù)測(cè)標(biāo)準(zhǔn)差，可以看出gz反演得到的gz分量的預(yù)測(cè)標(biāo)準(zhǔn)差小于其他組合反演得到的gz分量的預(yù)測(cè)標(biāo)準(zhǔn)差，而gz反演中梯度分量的預(yù)測(cè)標(biāo)準(zhǔn)差則大于其他分量組合反演得到的梯度分量的預(yù)測(cè)標(biāo)準(zhǔn)差，這是因?yàn)樵诜囱萦?jì)算中，參與反演過(guò)程的數(shù)據(jù)擬合的更好.所以，將多個(gè)分量聯(lián)合反演，數(shù)據(jù)擬合的結(jié)果更好，反演結(jié)果更加合理.

圖20d—20f展示了gzz反演結(jié)果.相比gz反演結(jié)果，gzz反演結(jié)果更加緊密.由圖20d可知，巖蓋的中心埋深大約是400 m,其中高密度部分的在東西向上的長(zhǎng)度大約是1400 m.這一結(jié)果與Thompson等(1928)中結(jié)果相符.由圖20e和圖20f還可以判斷出，巖丘由東北到西南方向的延長(zhǎng)，這一特征與該地區(qū)斷層的走向吻合(Coker et al.，2007).所以， gzz反演結(jié)果是合理的.

對(duì)gz|gzz進(jìn)行反演計(jì)算，結(jié)果如圖20g—20i所示.對(duì)圖20d和圖20g進(jìn)行比較，發(fā)現(xiàn)圖20g更加緊密.對(duì)表2中g(shù)z|gzz反演得到的各分量的預(yù)測(cè)標(biāo)準(zhǔn)差(第3列)和gzz反演得到的各分量的預(yù)測(cè)標(biāo)準(zhǔn)差(第2列)進(jìn)行比較，可以得出gz|gzz反演得到的梯度分量的預(yù)測(cè)標(biāo)準(zhǔn)差增加了，這說(shuō)明加入gz數(shù)據(jù)能夠識(shí)別出與gz中包含信息不符的噪聲，得到的結(jié)果更加合理.

圖21a—21c展示了gxz|gyz|gzz反演結(jié)果.圖21d—21f展示了gz|gxz|gyz|gzz結(jié)果.比較圖21a和圖21d發(fā)現(xiàn)，圖21d在圖21a的基礎(chǔ)上更進(jìn)一步優(yōu)化，高密度部分更加連續(xù).圖21g—21i展示了gxy|gxz|gyy|gyz|gzz反演結(jié)果，對(duì)比圖21g和圖21d,模型高密度部分更加緊密，更加連續(xù)，模型得到了進(jìn)一步的優(yōu)化.

圖21j—21l給出了gz|gxy|gxz|gyy|gyz|gzz反演結(jié)果，該結(jié)果與gxy|gxz|gyy|gyz|gzz反演結(jié)果類似.但是，對(duì)比表2中由兩組分量組合反演結(jié)果計(jì)算得到的各分量的預(yù)測(cè)標(biāo)準(zhǔn)差(第6行和第7行)，可以發(fā)現(xiàn)，gz|gxy|gxz|gyy|gyz|gzz反演得到的各分量的預(yù)測(cè)標(biāo)準(zhǔn)差更小，說(shuō)明gz|gxy|gxz|gyy|gyz|gzz反演中各數(shù)據(jù)擬合更好，反演結(jié)果更加合理.文中選擇gz|gxy|gxz|gyy|gyz|gzz反演結(jié)果作為最終預(yù)測(cè)的模型.

圖22給出了gz|gxy|gxz|gyy|gyz|gzz的橫切面，橫切面的深度從200 m變化到600 m，間隔為50 m.由圖中可以看出，模型的南部的深度大約為200 m，由南向北，模型深度增大，這與Thompson等(1928)中的鉆井?dāng)?shù)據(jù)相吻合，資料顯示Vinton巖丘南部深度大約為200 m，向北部深度達(dá)到315 m.除此之外，模型的底部大約是600 m.Oliveira等(2013)給出的模型的深度是460 m，Geng等(2014)給出的深度是650 m.如果想要進(jìn)一步確認(rèn)準(zhǔn)確的深度，需要更多的資料來(lái)約束反演，這里，我們認(rèn)為，得出的模型是合理的.

圖18　Vinton鹽丘的觀測(cè)數(shù)據(jù)(a) gxx， (b) gxy， (c) gxz， (d) gyy， (e) gyz， (f) gz， (g) gzz. (g)中的黑線表示圖19,圖20,圖21中垂直剖面的位置.Fig.18　Observed data of Vinton salt dome(a) gxx, (b) gxy, (c) gxz, (d) gyy, (e) gyz, (f)gz, (g) gzz. Black line in (g) shows the location of the vertical profile in Figs.19,20, and 21.

圖20　gz反演結(jié)果(a，b，c)，gzz反演結(jié)果(d，e，f)，gz|gzz反演結(jié)果(g，h，i)，其中左側(cè)一列為垂直剖面，剖面位置如圖18g中黑線所示，中間一列為橫切面，深度為z=250 m，右側(cè)一列為橫切面，深度為z=400 mFig.20　gz inversion result is shown in (a,b,c), gzz inversion result is shown in (d,e,f), gz|gzz inversion result is shown in (g,h,i). The Left column is vertical profile, the location of the profile is shown in Fig.18g. The central column is the depth section at z=250 m. The right column is the depth section at z=400 m

gxxgxygxzgyygyzgzzgzgz6.4004.5077.9667.7418.93711.8580.289gzz2.2451.5793.0803.2272.7923.3920.432gz|gzz2.3821.5983.2072.2992.8883.5440.430gxz|gyz|gzz2.2011.5632.9082.1082.7303.3140.459gz|gxz|gyz|gzz2.2391.5672.9072.1272.7233.3720.454gxy|gxz|gyy|gyz|gzz2.1161.4182.6501.7652.3652.9470.455gz|gxy|gxz|gyy|gyz|gzz2.1011.4052.6201.7362.3322.9100.450

6結(jié)論

本文首先提出了一種可以快速計(jì)算重力和重力梯度靈敏度矩陣的方法，對(duì)于地下規(guī)則劃分的情況，采用這一方法計(jì)算，能夠大大地縮短計(jì)算時(shí)間.針對(duì)各分量聯(lián)合反演數(shù)據(jù)量增大的問(wèn)題，文中選擇有限內(nèi)存BFGS擬牛頓法求解目標(biāo)函數(shù)，這一方法對(duì)內(nèi)存要求相對(duì)較小，能夠使目標(biāo)函數(shù)較快地收斂.為了克服趨膚效應(yīng)，文中介紹了加權(quán)方式和深度加權(quán)函數(shù)，并對(duì)目標(biāo)函數(shù)的梯度的表達(dá)式進(jìn)行了推導(dǎo).本文得出，由單一分量反演獲得的異常體埋深的估計(jì)值與異常體的實(shí)際埋深存在一定的偏差，但是將估計(jì)值應(yīng)用到實(shí)際反演計(jì)算中，得到的反演結(jié)果與實(shí)際模型是一致的.而且，當(dāng)估計(jì)值在一定范圍內(nèi)波動(dòng)時(shí)，得到的反演結(jié)果也是可信的.通過(guò)第一組組合模型試驗(yàn)，發(fā)現(xiàn)正則化參數(shù)取1是可行的，不同噪聲水平的數(shù)據(jù)對(duì)正則化參數(shù)的取值影響不大.通過(guò)第二組組合模型試驗(yàn)，發(fā)現(xiàn)將重力數(shù)據(jù)和重力梯度數(shù)據(jù)聯(lián)合反演能夠在一定程度上提高反演結(jié)果，但是結(jié)果的提升并不明顯.這可能是因?yàn)橹亓?shù)據(jù)包含較多的低頻信息，而重力梯度數(shù)據(jù)包含較多的高頻信息，因此，重力數(shù)據(jù)包含更多的深部信息，重力梯度數(shù)據(jù)包含更多的淺部信息.所以將重力數(shù)據(jù)和重力梯度數(shù)據(jù)聯(lián)合反演淺部地質(zhì)體，得到的反演結(jié)果與重力梯度數(shù)據(jù)反演結(jié)果相比，提升并不明顯.在后續(xù)的研究中，我們將針對(duì)深部地質(zhì)體進(jìn)行討論.將該方法應(yīng)用到美國(guó)路易斯安那州的Vinton巖丘的FTG數(shù)據(jù)和重力數(shù)據(jù)，發(fā)現(xiàn)通過(guò)聯(lián)合反演，模型得到了一定程度的提升，這說(shuō)明在實(shí)際數(shù)據(jù)應(yīng)用中，對(duì)于淺部地質(zhì)體，聯(lián)合反演是有一定意義的.反演得到的結(jié)果和已知的地質(zhì)資料吻合，證明該方法可以用于實(shí)際數(shù)據(jù)的處理.

圖19　gzz反演結(jié)果的垂直剖面，剖面位置如圖18g中黑線所示.由圖中可知，巖蓋位于地下200 m以下Fig.19　Vertical profile of gzz inversion result, the location of the profile is shown in Fig.18g. The cap-rock is located 200m below the surface

致謝感謝Bell Geospace公司提供數(shù)據(jù).

References

Avdeev D, Avdeeva A. 2009. 3D magnetotelluric inversion using a limited-memory quasi-Newton optimization.Geophysics, 74(3): F45-F57.

Avdeeva A, Avdeev D. 2006. A limited-memory quasi-Newton inversion for 1D magnetotellurics.Geophysics, 71(5): G191-G196.

Capriotti J, Li Y G. 2014. Gravity and gravity gradient data: Understanding their information content through joint inversions.∥ 2014 SEG Annual Meeting. Denver, Colorado, USA: SEG, 1329-1333.

Chasseriau P, Chouteau M. 2003. 3D gravity inversion using a model of parameter covariance.JournalofAppliedGeophysics, 52(1): 59-74.Coker M O, Bhattacharya J P, Marfurt K J. 2007. Fracture patterns within mudstones on the flanks of a salt dome: Syneresis or slumping? Gulf Coast Association of Geological Societies Transactions, 57: 125-137.

Commer M. 2011. Three-dimensional gravity modeling and focusing inversion using rectangular meshes.GeophysicalProspecting, 59(5): 966-979.

Commer M, Newman G A, William K H, et al. 2011. 3D induced-polarization data inversion for complex resistivity.Geophysics, 76(3): F157-F171.

Ennen C. 2012. Mapping gas-charged fault blocks around the Vinton Salt Dome, Louisiana using gravity gradiometry data [Master′s thesis]. Houston: University of Houston. Farquharson C G, Oldenburg D W. 2004. A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems.Geophysics, 156(3): 411-425. Forsberg R. 1984. A study of terrain corrections, density anomalies, and geophysical inversion methods in gravity field modeling.∥ Report 355, Department of Geodetic Science and Surveying, The Ohio Stata University.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. Haber E. 2005. Quasi-Newton methods for large-scale electromagnetic inverse problems.InverseProblems, 21(1): 305-323.Last B J, Kubik K. 1983. Compact gravity inversion.Geophysics, 48(6): 713-721.

Li X, Chouteau M. 1998. Three-dimensional gravity modeling in all space.SurveysinGeophysics, 19(4): 339-368.

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

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.∥ SEG Int′l Exposition and Annual Meeting. San Antonio, Texas.

Li Y G, Oldenburg D W. 2003. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method.GeophysicalJournalInternational, 152(2): 251-265.

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.

Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data.ChineseJ.Geophys. (in Chinese), 56(12): 4278-4287, doi: 10.6038/cjg20131230.

Liu Y P, Wang Z W, Du X J, et al. 2013. 3D constrained inversion of gravity data based on the Extrapolation Tikhonov regularization algorithm.ChineseJ.Geophys. (in Chinese), 56(5): 1650-1659, doi: 10.6038/cjg20130522. Newman G A, Boggs P T. 2004. Solution accelerators for large-scale three-dimensional electromagnetic inverse problems.InverseProblems, 20(6): S151-S170. Ni Q, Yuan Y. 1997. A subspace limited memory quasi-Newton algorithm for large-scale nonlinear bound constrained optimization.MathematicsofComputation, 66(220): 1509-1520. Nocedal J, Wright S J. 1999. Numerical Optimization. New York: Springer.

Oliveira V C Jr, Barbosa V C F. 2013. 3-D radial gravity gradient inversion.GeophysicalJournalInternational, 195(2): 883-902. Shamsipour P, Denis M, Michel C, et al. 2010. 3D stochastic inversion of gravity data using cokriging and cosimulation.Geophysics, 75(1): I1-I10.Plessix R E, Mulder W A. 2008. Resistivity imaging with controlled-source electromagnetic data: Depth and data weighting.InverseProblems, 24(3): 1-22. Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images.Geophysics, 64(3): 874-887.

Portniaguine O, Zhdanov M S. 2002. 3-D magnetic inversion with data compression and image focusing.Geophysics, 67(5): 1532-1541.

圖22　gz|gxy|gxz|gyy|gyz|gzz 反演模型的橫切面，切面深度變化范圍是200 m到600 m，間隔為50 mFig.22　Depth sections of gz|gxy|gxz|gyy|gyz|gzz inversion results. The depth ranges from 200 m to 600 m. The interval is 50 m

Thompson S A, Eichelberger O H. 1928. Vinton salt dome, Calcasieu Parish, Louisiana.AAPGBulletin, 12(4): 385-394. Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. Washington, D.C.: W. H. Winston and Sons.

Wu L, Ke X P, Hsu H, et al. 2013. Joint gravity and gravity gradient inversion for subsurface object detection.IEEEGeoscienceandRemoteSensingLetters, 10(4): 865-849.

Yao C L, Hao T Y, Guan Z N, et al. 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.

Zhang L L, Koyama T, Utada H, et al. 2012. A regularized three-dimensional magnetotelluric inversion with a minimum gradient support constraint.GeophysicalJournalInternational, 189(1): 296-316.

Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. New York: Elsevier.

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

Zhdanov M S. 2009. New advances in regularized inversion of gravity and electromagnetic data.GeophysicalProspecting, 57(4): 463-478.

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

劉云鶴, 殷長(zhǎng)春. 2013. 三維頻率域航空電磁反演研究. 地球物理學(xué)報(bào), 56(12): 4278-4287, doi: 10.6038/cjg20131230.

劉銀萍, 王祝文, 杜曉娟等. 2013. 基于Extrapolation Tikhonov正則化算法的重力數(shù)據(jù)三維約束反演. 地球物理學(xué)報(bào), 56(5): 1650-1659, doi: 10.6038/cjg20130522.

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

(本文編輯何燕)

基金項(xiàng)目國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)課題(2014AA06A613)資助.

作者簡(jiǎn)介秦朋波，1988年生，男，在讀博士，主要從事重磁和梯度數(shù)據(jù)反演研究. E-mail:qin-pengbo@163.com

doi:10.6038/cjg20160624 中圖分類號(hào)P631

收稿日期2015-10-09，2016-02-18收修定稿

Integrated gravity and gravity gradient data focusing inversion

QIN Peng-Bo, HUANG Da-Nian

CollegeofGeo-explorationSciencesandTechnology,JilinUniversity,Changchun130000,China

AbstractGravity data contains a wealth of low frequency contents and gravity gradient data contains many high frequency contents. Thus, in theory, we can get a more reliable result through integrated gravity and gravity gradient inversion. Here, we propose a method to integrate gravity and gravity gradient data in inversion. The computation time and requirement for memory of the inversion are increased with multiple components included. Thus, we present a method to calculate sensitivity matrixes of different components to reduce the computational time. We adopt limited-memory BFGS quasi-Newton algorithm to solve the inverse problem. It uses curvature information from only the most recent iterations to construct the Hessian approximation. The requirement for storage is reduced in this way. A weighting scheme for resolution enhancement at depth is introduced through the re-weighted method. We estimate the depth of the anomalous body by single component inversion. Then the depth information is incorporated into the depth weighting functional. At last, we adopt re-weighted method to combine the depth weighting functional with the objective functional. We use a synthetic example to demonstrate the process. Although the estimate of depth is not accurate, we can get a more accurate inversion result by combine the depth information into the depth weighting function and apply it in inversion. We adopt a synthetic example to explore the advantages of integrated gravity and gravity gradient inversion. The results show that integrated gravity data and gravity gradient data in inversion, the addition noise which is inconsistent with the corresponding components can be identified. The recovered model is more reasonable. However, for different component combinations, the inversion results are similar, which indicates that the improvement of the recovered model is small. At last, we apply the method to real data of the Vinton salt dome, Louisiana, USA. The results indicate that the recovered model is improved through integrated gravity and gravity gradient inversion.

KeywordsGravity and gravity gradient forward; Integrated inversion of gravity and gravity gradient; Limited-memory BFGS quasi Newton algorithm; Depth weighting functional; Minimum gradient support functional

秦朋波， 黃大年. 2016. 重力和重力梯度數(shù)據(jù)聯(lián)合聚焦反演方法. 地球物理學(xué)報(bào)，59(6):2203-2224,doi:10.6038/cjg20160624.Qin P B, Huang D N. 2016. Integrated gravity and gravity gradient data focusing inversion.ChineseJ.Geophys. (in Chinese),59(6):2203-2224,doi:10.6038/cjg20160624.