黃天統(tǒng)，彭新發(fā)，朱自強(qiáng)

(核工業(yè)二三〇研究所，湖南 長(zhǎng)沙 410007)

0 引言

作為一種高精度的地球物理勘探方法，在實(shí)際勘探過(guò)程中重力梯度張量測(cè)量的是重力位的二階導(dǎo)數(shù)。與傳統(tǒng)的重力勘探相比，重力梯度張量測(cè)量對(duì)地下介質(zhì)密度變化反映更加靈敏，具有更好的抗干擾能力，能夠更加直接地反映地下密度異常體的水平邊界位置[1-3]。在國(guó)外，重力梯度張量勘探技術(shù)已經(jīng)較為廣泛地應(yīng)用于海洋領(lǐng)域、航空領(lǐng)域以及衛(wèi)星重力領(lǐng)域。

結(jié)構(gòu)化網(wǎng)格是指網(wǎng)格區(qū)域內(nèi)所有的內(nèi)部點(diǎn)都具有相同的毗鄰單元。由于非結(jié)構(gòu)網(wǎng)格具有生成速度快、網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單、網(wǎng)格的質(zhì)量好等優(yōu)點(diǎn)[4]，廣泛地應(yīng)用于形狀規(guī)則的物體剖分。但是隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，結(jié)構(gòu)化網(wǎng)格的弊端越來(lái)越明顯，比如：待剖分物體的幾何形狀越來(lái)越復(fù)雜，結(jié)構(gòu)化網(wǎng)格剖分帶來(lái)了的剖分誤差也越來(lái)越大，并且節(jié)點(diǎn)和網(wǎng)格單元在高質(zhì)量情況下不能合理地呈梯度狀分布在剖分的網(wǎng)格中[5]。

與結(jié)構(gòu)化網(wǎng)格相比，非結(jié)構(gòu)化網(wǎng)格彌補(bǔ)了其不能解決所有形狀和任意連通區(qū)域網(wǎng)格剖分的不足；對(duì)于網(wǎng)格單元和節(jié)點(diǎn)的編號(hào)，非結(jié)構(gòu)化網(wǎng)格并沒(méi)有固定的規(guī)則，兩個(gè)相鄰單元格編號(hào)對(duì)應(yīng)的網(wǎng)格單元不一定相鄰，這樣使得非結(jié)構(gòu)化網(wǎng)格能夠十分方便的控制單元和節(jié)點(diǎn)的分布，具有很強(qiáng)的適應(yīng)性，能夠更好地剖分復(fù)雜模型。

在地球物理數(shù)值模擬過(guò)程中，高效地?cái)M合場(chǎng)源模型及地形分布將更加有利于模型的計(jì)算，非結(jié)構(gòu)化網(wǎng)格與生俱來(lái)的優(yōu)勢(shì)使其廣泛地應(yīng)用于地球物理數(shù)值模擬中。湯井田[6]、任政勇[7]在基于非結(jié)構(gòu)化網(wǎng)格的情況下實(shí)現(xiàn)了2.5維的自適應(yīng)有限元直流電阻率的正演模擬；孫麗影[8]、劉長(zhǎng)生[9]、童孝忠[10]通過(guò)非結(jié)構(gòu)化網(wǎng)格進(jìn)行局部加密完成了大地電磁二維、三維有限元數(shù)據(jù)模擬；趙曉博利用Delaunay三角化實(shí)現(xiàn)了中心回線瞬變電磁2.5維的正演模擬；曾思紅[11]、曹書(shū)錦[12]利用有限單元法在非結(jié)構(gòu)化網(wǎng)格下對(duì)重力梯度張量進(jìn)行了二維、三維正演，Jahandari[13]則是運(yùn)用有限體積法實(shí)現(xiàn)了基于非結(jié)構(gòu)化網(wǎng)格的二維、三維的重力高精度正演模擬，取得了較好的效果；Okabe[14]推導(dǎo)了任意多邊形和多面體網(wǎng)格下的重磁一階、二階導(dǎo)數(shù)的二維、三維解析表達(dá)式，具有很高的計(jì)算精度。

目前，對(duì)于重力梯度張量數(shù)據(jù)反演，大多都是基于長(zhǎng)方形或者長(zhǎng)方體等結(jié)構(gòu)化網(wǎng)格進(jìn)行的，但由于結(jié)構(gòu)化網(wǎng)格的幾何適應(yīng)能力差，在較復(fù)雜地質(zhì)體或地形的情況下會(huì)帶來(lái)較大的剖分誤差，從而使得反演效果有所下降。與結(jié)構(gòu)化網(wǎng)格相比，非結(jié)構(gòu)化網(wǎng)格能夠減少模型的剖分誤差，更好地?cái)M合地下密度異常體的分布和地表地形的變化，因此，研究基于非結(jié)構(gòu)化網(wǎng)格的重力梯度張量反演算法，可以提高重力梯度張量反演過(guò)程中的計(jì)算精度?，F(xiàn)在，根據(jù)非結(jié)構(gòu)化網(wǎng)格進(jìn)行的地球物理反演主要還是基于二維三角形網(wǎng)格，elièvre[15]、吉日嘎拉圖[16]利用地震初至波完成了二維慢度的反演，并完成了非結(jié)構(gòu)化網(wǎng)格重震二維聯(lián)合反演，使得反演的縱橫向分辨率都很高，降低了反演結(jié)果的多解性問(wèn)題。對(duì)于位場(chǎng)三維非結(jié)構(gòu)化網(wǎng)格的反演，國(guó)內(nèi)外文獻(xiàn)中很少有涉及到。

1 正反演算法

1.1 基于非結(jié)構(gòu)化網(wǎng)格重力梯度張量正演算法

重力梯度張量是重力位的二階導(dǎo)數(shù)，在笛卡爾坐標(biāo)系下，根據(jù)牛頓萬(wàn)有引力定律，一個(gè)剩余密度為ρ，觀測(cè)點(diǎn)為P(p0,q0,r0)，位于Q(p,q,r)地下異常體ν在地下產(chǎn)生的重力位V有：

(1)

其中：G為萬(wàn)有引力常量，u=-(x2+y2+z2)-1/2,x=p-p0、y=q-q0、z=r-r0。則重力梯度張量表達(dá)式寫(xiě)成矩陣形式：

(2)

其中：Γ表示全重力梯度張量。重力梯度張量共有9個(gè)分量，由于在無(wú)源空間里重力的旋度和散度都為0，有Vxx+Vyy+Vzz=0，Vyx=Vxy，Vzx=Vxz，Vzy=Vyz，因此重力梯度張量中只有5個(gè)獨(dú)立分量Vxx、Vxy、Vxz、Vyz、Vzz。

重力位在k方向上的一階導(dǎo)數(shù)：

(3)

其中：k是k方向法向量：kT=(kx,ky,kz)或者kT=(cos (x,k),cos (y,k),cos (z,k))。

將式(1)帶入式(3)得：

(4)

引入另一個(gè)方向向量l：lT=(lx,ly,lz)，對(duì)于不同的重力梯度張量分量，k、l的方向余弦取值如表1所示。對(duì)Vk求在l方向求導(dǎo)，得到重力位的二階導(dǎo)數(shù)：

(5)

表1 方向余弦的取值與重力梯度張量各分量的關(guān)系

對(duì)式(5)利用散度定理得：

(6)

其中：S為地下異常體ν的表面，n為面S的法向量，指向異常體外部。

在四面體網(wǎng)格的形式下,式(6)的積分可以寫(xiě)成每個(gè)網(wǎng)格的積分總和的形式：

(7)

其中：

(8)

為了快速地求解出面Si上Ki的值，采用了坐標(biāo)選擇的方法。旋轉(zhuǎn)分為兩步：首先，以z軸為中心逆時(shí)針[15]旋轉(zhuǎn)xOy平面(笛卡爾坐標(biāo)系)，直到x方向與面的外法向量在xOy平面上的投影平行，旋轉(zhuǎn)的角度為θ，得到最終的坐標(biāo)系統(tǒng)(X,Y,Z)。然后，以Y軸為中心逆時(shí)針選擇x′Oz軸，直到z方向與面的外法向量平行，旋轉(zhuǎn)的角度為φ，得到最終的坐標(biāo)系統(tǒng)(X,Y,Z)，如圖1所示。

圖1 面坐標(biāo)系統(tǒng)旋轉(zhuǎn)示意

坐標(biāo)的旋轉(zhuǎn)轉(zhuǎn)換可以寫(xiě)成：

(9)

其中：0≤θ<2π,0≤φ≤π。通過(guò)旋轉(zhuǎn)，面Si上的3個(gè)頂點(diǎn)的Z值相等，面Si的外法向量ni可以寫(xiě)成：

(10)

在(X,Y,Z)坐標(biāo)系統(tǒng)下，式(8)可以寫(xiě)成：

(11)

其中：

(12)

通過(guò)對(duì)Si面上的每條邊進(jìn)行二維坐標(biāo)旋轉(zhuǎn)，可以將式(11)的面積分轉(zhuǎn)換為線積分。以原點(diǎn)為中心，逆時(shí)針旋轉(zhuǎn)X軸、Y軸，直到X方向與Si面上的j邊重合，得到新的坐標(biāo)系統(tǒng)(ξ,η)旋轉(zhuǎn)的角度為ψ，0≤ψ<2π，如圖2所示。坐標(biāo)的旋轉(zhuǎn)轉(zhuǎn)換可以寫(xiě)成：

(13)

通過(guò)旋轉(zhuǎn)，Si面上j邊的兩個(gè)端點(diǎn)的η值相等，j邊的法向量sj為：

圖2 邊坐標(biāo)系統(tǒng)旋轉(zhuǎn)示意

(14)

那么，式(11)就可以寫(xiě)成：

(15)

從而得到：

(16)

對(duì)式(16)求積分得到：

ln[ξ+(ξ2+η2+Z2)1/2]·

(17)

最終得到重力位二階導(dǎo)數(shù)的表達(dá)式:

(18)

如果Z或者cosψ的值為0，則式(17)的后面一項(xiàng)為0，但是此時(shí)式(17)的第一項(xiàng)的值不一定一直存在，例如η和Z的值都趨近于0，且ξ不為正數(shù)。當(dāng)ξj和ξj+1都為負(fù)時(shí)，通過(guò)式(19)可以避免，對(duì)于其他情形可用一個(gè)較小的值代替。

(19)

1.2 反演的基本理論

由于重力梯度張量反問(wèn)題是一個(gè)不適定問(wèn)題，為了得到一個(gè)穩(wěn)定的可行解，李耀國(guó)等人[16]在吉洪諾夫正則化理論的基礎(chǔ)上提出了一個(gè)能夠廣泛應(yīng)用于地球物理反問(wèn)題的廣義目標(biāo)函數(shù)，其中數(shù)據(jù)擬合函數(shù)為:

其中：Wd=diag{1/ε1,…,1/εN}，εi為第i個(gè)觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)差。通常假設(shè)實(shí)際數(shù)據(jù)的誤差是相互獨(dú)立且滿足高斯分布的，實(shí)測(cè)得到重力梯度張量數(shù)據(jù)dobs=F(m)+δ，δ表示的是觀測(cè)點(diǎn)上的隨機(jī)噪聲。由于隨機(jī)噪聲δ是滿足高斯分布N(0,ε2)的，ε為觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)差。根據(jù)最小二乘擬合，那么反演最終的數(shù)據(jù)擬合差為：

(21)

其中：N為測(cè)點(diǎn)數(shù)。

在最小化目標(biāo)函數(shù)的過(guò)程中，在目標(biāo)函數(shù)中引入的模型目標(biāo)函數(shù)為：

(22)

φm(m)=αs‖Ws(m-mref)‖2+

αx‖Wx(m-mref)‖2+αy‖Wy(m-mref)‖2+

αz‖Wz(m-mref)‖2

(23)

其中：αs、αx、αy、αz為每項(xiàng)的加權(quán)系數(shù)，Ws為粗糙度矩陣，Wx、Wy、Wz分別為x、y、z三個(gè)方向上的差分算子。令：

(24)

然而，在非結(jié)構(gòu)化網(wǎng)格下，差分矩陣Wx、Wy、Wz的求解不能像在結(jié)構(gòu)化網(wǎng)格情況下那樣簡(jiǎn)單的求解。為了解決差分矩陣的計(jì)算問(wèn)題，在Lelièvre[15]提出的二維求解Wm的計(jì)算方法基礎(chǔ)上進(jìn)行衍生，得到三維情況下兩個(gè)相鄰網(wǎng)格Wm的計(jì)算方法。

(25)

其中：V1、V2為兩個(gè)相鄰四面體網(wǎng)格的體積，R為兩個(gè)四面體網(wǎng)格的中心距離。

結(jié)合數(shù)據(jù)擬合函數(shù)式(1)和模型目標(biāo)函數(shù)式(4)可以得到廣義的反演目標(biāo)函數(shù)：

1.3 深度加權(quán)函數(shù)

為了使得反演結(jié)果能夠反演出深度上的結(jié)構(gòu)，提高深度方向上的分辨率，李耀國(guó)等人[19]引入了一個(gè)深度加權(quán)函數(shù)來(lái)減小核矩陣在深度方向上的衰減：

(27)

其中：z是網(wǎng)格單元的中心埋深，z0是個(gè)常數(shù)，β是一個(gè)接近于3的常數(shù)，在已知網(wǎng)格剖分方式和觀測(cè)面的高度后，通過(guò)調(diào)節(jié)z0和β的值，來(lái)盡可能地減少核函數(shù)在深度方向上的衰減，從而使得反演結(jié)果能夠準(zhǔn)確地分布在對(duì)應(yīng)的深度上。式(24)結(jié)合深度加權(quán)函數(shù)得到

(28)

將深度加權(quán)函數(shù)結(jié)合式(22)得到帶深度加權(quán)的模型目標(biāo)函數(shù)：

φm(m)=‖Wm(m-mref)‖2

(29)

1.4 多分量聯(lián)合反演

多分量的聯(lián)合反演的目標(biāo)函數(shù)只需要對(duì)式(26)中的數(shù)據(jù)擬合項(xiàng)進(jìn)行適當(dāng)?shù)男薷模磳⒏鱾€(gè)分量的數(shù)據(jù)擬合項(xiàng)結(jié)合起來(lái)形成一個(gè)總的數(shù)據(jù)擬合項(xiàng)。結(jié)合模型目標(biāo)函數(shù)，對(duì)于重力梯度張量多分量聯(lián)合反演的目標(biāo)函數(shù)可以寫(xiě)成：

a‖Wm(m-mref)‖2，

(30)

α,β=x,y,z

2 算法驗(yàn)算

圖3所示為單個(gè)長(zhǎng)方體模型，場(chǎng)源空間大小為1 000 m×1 000 m×500 m，其內(nèi)埋藏著一個(gè)大小400 m×400 m×200 m的存在密度異常的長(zhǎng)方體，長(zhǎng)方體頂部埋深為150 m，底部埋深為350 m，水平中心位置在地表的投影坐標(biāo)為(500 m,500 m)，長(zhǎng)方體剩余密度為1.0 g/cm3，觀測(cè)位置位于地表，x、y方向測(cè)線的長(zhǎng)度1 000 m、點(diǎn)距50 m。根據(jù)已剖分的網(wǎng)格單元，應(yīng)用前文推導(dǎo)的公式計(jì)算出來(lái)的重力梯度張量各分量的正演結(jié)果如圖4所示。

圖3 單個(gè)長(zhǎng)方體模型非結(jié)構(gòu)化網(wǎng)格剖分示意

圖4 單個(gè)長(zhǎng)方體模型重力梯度張量正演結(jié)果

從圖4可以看出，重力梯度張量各個(gè)分量的正演結(jié)果都能夠明顯地反映地下密度異常體的邊界位置。Vxy的4個(gè)極值點(diǎn)對(duì)應(yīng)了長(zhǎng)方體4個(gè)頂點(diǎn)的水平投影位置；Vxx、Vyy的極小值位置對(duì)應(yīng)的是長(zhǎng)方體水平中心位置，并且等值線變化最密集的部分正好是長(zhǎng)方體的4條邊界位置；Vxz、Vyz的極值位置對(duì)應(yīng)的是長(zhǎng)方體邊界的水平中心位置；Vzz的等值線圖類(lèi)似一系列的同心矩形，等值線在地下異常體密度變化邊界分布密集，能夠很好地確定密度變化，極大值的中心位置對(duì)應(yīng)的是長(zhǎng)方體的中心位置。

將單個(gè)長(zhǎng)方體模型計(jì)算出來(lái)的重力梯度張量正演結(jié)果，利用正則化方法并結(jié)合深度加權(quán)函數(shù)對(duì)重力梯度張量各分量進(jìn)行反演，其中初始參考模型mref=0，深度加權(quán)矩陣中的β取值為 3。各分量反演結(jié)果如圖5所示。

從圖5可以看出，反演結(jié)果對(duì)異常頂部位置的分辨率比底部分辨率要高，這是由于位場(chǎng)本身的特征決定的。對(duì)比各個(gè)分量之間反演效果，各分量反演出來(lái)的剩余密度值與真實(shí)值比較接近，但是都出現(xiàn)了一定的負(fù)值，其中，Vxz、Vyy、Vyz反演出來(lái)的負(fù)值較大；各分量的反演結(jié)果在深度方向上延伸較大；Vzz反演出的剩余密度更接近真實(shí)值，反演的剩余密度極大值位置與真實(shí)位置非常吻合,Vxx、Vxy、Vxz、Vyy、Vyz、Vzz反演的密度值與理論模型密度的均方誤差分別是0.214 8、0.205 3、0.201 5、0.217 2、0.213 1、0.192 7，說(shuō)明了對(duì)于單分量的反演，基于非結(jié)構(gòu)化網(wǎng)格的重力梯度張量反演方法是十分有效的，算法是可行的。

圖5 單個(gè)長(zhǎng)方體模型單分量反演結(jié)果

為了充分利用重力梯度張量各分量包含的密度、空間信息，對(duì)重力梯度張量的獨(dú)立分量進(jìn)行多分量的聯(lián)合反演，結(jié)果如圖6所示。對(duì)比圖3與圖6可知，多分量聯(lián)合反演出來(lái)的密度值與理論模型的密度的均方誤差為0.190 9，比任一單分量反演的均方誤差都小，且利用多分量進(jìn)行聯(lián)合反演，其結(jié)果與理論模型極其吻合，這也證明了多分量聯(lián)合反演的正確性。與單分量反演相比，聯(lián)合反演出來(lái)的剩余密度值更加接近真實(shí)值，反演結(jié)果在水平方向和深度方向上都更加集中收斂于異常體的中心位置，利用多個(gè)獨(dú)立分量進(jìn)行聯(lián)合反演能夠更加準(zhǔn)確反演出地下介質(zhì)的剩余密度及埋藏位置，再次驗(yàn)證了算法的正確性。

圖6 單個(gè)長(zhǎng)方體模型多分量聯(lián)合反演結(jié)果

3 Y型巖脈模型

場(chǎng)源空間大小為1 000 m×1 000 m×500 m，如圖7所示，場(chǎng)源空間存在一個(gè)Y型巖脈，其剩余密度都為1.0 g/cm3，左側(cè)板狀體頂部埋深為100 m，頂部中心在地表的投影位置為(200 m,500 m)，頂部尺寸為200 m×400 m，底部埋深為300 m，底部中心在地表的投影位置為(300 m,500 m)，頂部尺寸為200 m×400 m。右側(cè)板狀體頂部埋深為100 m，頂部中心在地表的投影位置為(700 m,500 m)，頂部尺寸為200 m×400 m，底部埋深為400 m，底部中心在地表的投影位置為(500 m,500 m)，頂部尺寸為200 m×400 m。x、y方向點(diǎn)距為50 m，測(cè)線長(zhǎng)度為1 000 m。利用非結(jié)構(gòu)化網(wǎng)格剖分軟件將場(chǎng)源空間剖分，網(wǎng)格數(shù)目為9 408。

加入5%的高斯噪聲后，利用重力梯度張量的5個(gè)相互獨(dú)立分量進(jìn)行多分量的聯(lián)合反演，多分量聯(lián)合反演的密度值與理論模型密度的均方誤差為0.2347，反演結(jié)果如圖8所示。

對(duì)比圖7與圖8可知，利用多分量聯(lián)合反演的結(jié)果與理論模型相似，反演的密度異常體能夠很好地確定異常的位置、傾向。

4 結(jié)論

1)筆者利用吉洪諾夫正則化方法結(jié)合李耀國(guó)提出的廣義目標(biāo)函數(shù)將重力梯度張量各個(gè)分量正演結(jié)果分別進(jìn)行了單分量反演；其次，為了充分利用重力梯度張量各分量數(shù)據(jù)使得反演結(jié)果更加準(zhǔn)確，將5個(gè)相互獨(dú)立的重力梯度張量分量進(jìn)行了聯(lián)合反演；最后，用長(zhǎng)方體模型驗(yàn)證了算法的正確性，與長(zhǎng)方體網(wǎng)格下的反演結(jié)果作對(duì)比說(shuō)明非結(jié)構(gòu)化網(wǎng)格下反演的優(yōu)勢(shì)，并通過(guò)Y型巖脈模型來(lái)說(shuō)明方法的適應(yīng)性。得出以下結(jié)論： 1) 對(duì)重力梯度張量各個(gè)分量的正演結(jié)果進(jìn)行吉洪諾夫正則化反演，反演過(guò)程引入了深度加權(quán)函數(shù)，反演結(jié)果能夠較好地確定地下異常體的水平、深度位置，反演剩余密度值也與理論值比較接近，驗(yàn)證了方法的正確性。

圖7 Y型巖脈模型非結(jié)構(gòu)化網(wǎng)格剖分示意

圖8 Y型巖脈模型多分量聯(lián)合反演結(jié)果

2) 與單分量的重力梯度張量反演相比，多分量的聯(lián)合反演更加充分地利用了各個(gè)分量的信息，使得反演結(jié)果能夠更加吻合初始模型。

3) 通過(guò)Y型模型說(shuō)明了基于非結(jié)構(gòu)化網(wǎng)格重力梯度張量反演算法具有較強(qiáng)的適應(yīng)能力，能夠較好地反演復(fù)雜模型。