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

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

0 引言

作為一種高精度的地球物理勘探方法，在實際勘探過程中重力梯度張量測量的是重力位的二階導(dǎo)數(shù)。與傳統(tǒng)的重力勘探相比，重力梯度張量測量對地下介質(zhì)密度變化反映更加靈敏，具有更好的抗干擾能力，能夠更加直接地反映地下密度異常體的水平邊界位置[1-3]。在國外，重力梯度張量勘探技術(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)部點都具有相同的毗鄰單元。由于非結(jié)構(gòu)網(wǎng)格具有生成速度快、網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)簡單、網(wǎng)格的質(zhì)量好等優(yōu)點[4]，廣泛地應(yīng)用于形狀規(guī)則的物體剖分。但是隨著計算機(jī)技術(shù)的飛速發(fā)展，結(jié)構(gòu)化網(wǎng)格的弊端越來越明顯，比如：待剖分物體的幾何形狀越來越復(fù)雜，結(jié)構(gòu)化網(wǎng)格剖分帶來了的剖分誤差也越來越大，并且節(jié)點和網(wǎng)格單元在高質(zhì)量情況下不能合理地呈梯度狀分布在剖分的網(wǎng)格中[5]。

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

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

目前，對于重力梯度張量數(shù)據(jù)反演，大多都是基于長方形或者長方體等結(jié)構(gòu)化網(wǎng)格進(jìn)行的，但由于結(jié)構(gòu)化網(wǎng)格的幾何適應(yīng)能力差，在較復(fù)雜地質(zhì)體或地形的情況下會帶來較大的剖分誤差，從而使得反演效果有所下降。與結(jié)構(gòu)化網(wǎng)格相比，非結(jié)構(gòu)化網(wǎng)格能夠減少模型的剖分誤差，更好地擬合地下密度異常體的分布和地表地形的變化，因此，研究基于非結(jié)構(gòu)化網(wǎng)格的重力梯度張量反演算法，可以提高重力梯度張量反演過程中的計算精度?，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é)果的多解性問題。對于位場三維非結(jié)構(gòu)化網(wǎng)格的反演，國內(nèi)外文獻(xiàn)中很少有涉及到。

1 正反演算法

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

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

(1)

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

(2)

其中：Γ表示全重力梯度張量。重力梯度張量共有9個分量，由于在無源空間里重力的旋度和散度都為0，有Vxx+Vyy+Vzz=0，Vyx=Vxy，Vzx=Vxz，Vzy=Vyz，因此重力梯度張量中只有5個獨立分量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)

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

(5)

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

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

(6)

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

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

(7)

其中：

(8)

為了快速地求解出面Si上Ki的值，采用了坐標(biāo)選擇的方法。旋轉(zhuǎn)分為兩步：首先，以z軸為中心逆時針[15]旋轉(zhuǎn)xOy平面(笛卡爾坐標(biāo)系)，直到x方向與面的外法向量在xOy平面上的投影平行，旋轉(zhuǎn)的角度為θ，得到最終的坐標(biāo)系統(tǒng)(X,Y,Z)。然后，以Y軸為中心逆時針選擇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)換可以寫成：

(9)

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

(10)

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

(11)

其中：

(12)

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

(13)

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

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

(14)

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

(15)

從而得到：

(16)

對式(16)求積分得到：

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

(17)

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

(18)

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

(19)

1.2 反演的基本理論

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

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

(21)

其中：N為測點數(shù)。

在最小化目標(biāo)函數(shù)的過程中，在目標(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為每項的加權(quán)系數(shù)，Ws為粗糙度矩陣，Wx、Wy、Wz分別為x、y、z三個方向上的差分算子。令：

(24)

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

(25)

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

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

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

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

(27)

其中：z是網(wǎng)格單元的中心埋深，z0是個常數(shù)，β是一個接近于3的常數(shù)，在已知網(wǎng)格剖分方式和觀測面的高度后，通過調(diào)節(jié)z0和β的值，來盡可能地減少核函數(shù)在深度方向上的衰減，從而使得反演結(jié)果能夠準(zhǔn)確地分布在對應(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ù)只需要對式(26)中的數(shù)據(jù)擬合項進(jìn)行適當(dāng)?shù)男薷?，即將各個分量的數(shù)據(jù)擬合項結(jié)合起來形成一個總的數(shù)據(jù)擬合項。結(jié)合模型目標(biāo)函數(shù)，對于重力梯度張量多分量聯(lián)合反演的目標(biāo)函數(shù)可以寫成：

a‖Wm(m-mref)‖2，

(30)

α,β=x,y,z

2 算法驗算

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

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

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

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

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

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

圖5 單個長方體模型單分量反演結(jié)果

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

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

3 Y型巖脈模型

場源空間大小為1 000 m×1 000 m×500 m，如圖7所示，場源空間存在一個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方向點距為50 m，測線長度為1 000 m。利用非結(jié)構(gòu)化網(wǎng)格剖分軟件將場源空間剖分，網(wǎng)格數(shù)目為9 408。

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

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

4 結(jié)論

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

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

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

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

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