黃天統(tǒng),彭新發(fā),朱自強(qiáng)
(核工業(yè)二三〇研究所,湖南 長沙 410007)
作為一種高精度的地球物理勘探方法,在實際勘探過程中重力梯度張量測量的是重力位的二階導(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)中很少有涉及到。
重力梯度張量是重力位的二階導(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)
由于重力梯度張量反問題是一個不適定問題,為了得到一個穩(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ù):
為了使得反演結(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)
多分量的聯(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
圖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é)果
場源空間大小為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é)果與理論模型相似,反演的密度異常體能夠很好地確定異常的位置、傾向。
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ù)雜模型。