起伏地形下基于非結(jié)構(gòu)網(wǎng)格模型降維的重力三維密度反演

高新宇, 魯寶亮,2,3*, 安國(guó)強(qiáng), 巨鵬,2,3

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院, 西安 710054 2 海洋油氣勘探國(guó)家工程研究中心, 北京 100028 3 長(zhǎng)安大學(xué)西部礦產(chǎn)資源與地質(zhì)工程教育部重點(diǎn)實(shí)驗(yàn)室, 西安 710054

0 引言

重力三維密度反演是揭示地球內(nèi)部物質(zhì)結(jié)構(gòu)的重要地球物理方法.地下空間的剖分建模對(duì)反演結(jié)果具有重要的影響.直立長(zhǎng)方體剖分在勘探問(wèn)題中被廣泛采用(Nagy,1966; 姚長(zhǎng)利等,2003),但其在復(fù)雜地形條件下,近地表擬合存在較大誤差,進(jìn)而影響到正反演精度.因此,為了可以有效地逼近地形以及地下復(fù)雜巖體邊界,四面體非結(jié)構(gòu)化網(wǎng)格的重磁正反演受到廣泛關(guān)注(曹書錦等,2010; Jahandari and Farquharson,2013; Ren et al.,2017; Ma et al.,2022).

最初對(duì)于三維非結(jié)構(gòu)化網(wǎng)格,大多是用棱柱體進(jìn)行計(jì)算,對(duì)于正演問(wèn)題,Bhattacharyya(1964)提出了由矩形棱柱引起的重磁異常解析解,然后Hjelt(1972,1974)分別計(jì)算了均質(zhì)傾斜棱柱體的磁場(chǎng)表達(dá)式和重力異常的解析表達(dá)式.四面體或多面體的重力異常公式從棱柱體的公式發(fā)展而來(lái).Okabe(1979)利用散度定理,通過(guò)沿多面體模型邊緣的線積分計(jì)算,推導(dǎo)了均質(zhì)多面體在任意方向上重力異常的完整解析表達(dá)式,并分析了可能存在的奇點(diǎn)問(wèn)題.Pohánka(1988)也采用線積分法從多面體計(jì)算重力場(chǎng),描述了數(shù)值計(jì)算中遇到的問(wèn)題,并提出了解決方法.Singh和Guptasarma(2001)在均勻磁化多面體的計(jì)算方法上,提出了一種計(jì)算任意均勻密度多面體重力場(chǎng)的新方法,可以對(duì)所有的觀測(cè)點(diǎn)進(jìn)行.

三維密度反演是常用的重力資料定量解釋方法,但存在多解性和不穩(wěn)定等問(wèn)題,正則化方法是較為完善的解決方法之一,其中最常見的正則化方法是Tikhonov正則化方法(Li and Oldenburg,1998, 2003; 秦朋波和黃大年,2016; 高秀鶴和黃大年,2017; 李澤林等,2019; 王泰涵等,2020; Ma et al.,2022).對(duì)于復(fù)雜情況,反演必須能夠結(jié)合不同類型的先驗(yàn)知識(shí)和添加的約束條件,以便重建合理的反演模型.前人提出了多種約束條件,主要包括最小構(gòu)造約束(Last and Kubik,1983)、深度加權(quán)約束(Li and Oldenburg,1998,2003)、參考模型約束(陳石等,2014)、聚焦約束(Zhdanov et al.,2004)等.在高精度的地球物理反演中,由于地下剖分單元過(guò)多而使得矩陣過(guò)大,造成內(nèi)存不足,運(yùn)算時(shí)間長(zhǎng)等問(wèn)題,通過(guò)減小模型數(shù)量以及并行計(jì)算可以較為有效地解決這一問(wèn)題.例如通過(guò)邊界識(shí)別方法選取重點(diǎn)反演區(qū)域從而進(jìn)行模型降維以減小反演的復(fù)雜度(Yang et al.,2019).

在前人的研究基礎(chǔ)上,針對(duì)起伏地形及復(fù)雜地質(zhì)體重力三維密度反演問(wèn)題,本文提出采用四面體非結(jié)構(gòu)網(wǎng)格剖分地下空間,并利用垂向?qū)?shù)和相關(guān)系數(shù)進(jìn)行模型降維從而減小解的空間,聚焦目標(biāo)區(qū)域密度反演; 引入深度加權(quán)約束和物性范圍約束,構(gòu)建了基于模型降維的吉洪諾夫正則化目標(biāo)函數(shù),采用預(yù)條件共軛梯度法求解實(shí)現(xiàn)了重力三維密度反演.其中,為了提高反演效率在計(jì)算過(guò)程引入OpenMP并行算法.并研究了各項(xiàng)約束以及模型降維在四面體非結(jié)構(gòu)網(wǎng)格正則化反演中的效果.通過(guò)模型試驗(yàn)以及美國(guó)密蘇里州東南部氧化鐵礦床反演實(shí)踐,驗(yàn)證了方法的正確性和有效性.

1 正演方法

在直角坐標(biāo)系中,重力異常Δg的形式為

(1)

其中,ρ是場(chǎng)源密度,(x,y,z)為觀測(cè)點(diǎn)坐標(biāo),(ξ,η,ζ)為場(chǎng)源坐標(biāo).

對(duì)一個(gè)任意的四面體,Okabe(1979)根據(jù)散度定理給出了以三角面逼近任意多面體的重力異常計(jì)算方法,其基本思想是將多面體的體積分轉(zhuǎn)化為面積分,再轉(zhuǎn)化為線積分,從而完成重力異常的計(jì)算.對(duì)于一個(gè)任意的四面體(圖1),即分別計(jì)算四個(gè)三角面的面積分,將面積分又轉(zhuǎn)化為線積分,則完成了一個(gè)四面體的體積分計(jì)算.

圖1 四面體示意圖

對(duì)于將任意形體用非結(jié)構(gòu)網(wǎng)格四面體進(jìn)行剖分,計(jì)算所有四面體在場(chǎng)點(diǎn)產(chǎn)生的重力值然后求和,即得到任意形體的重力異常.假設(shè)將復(fù)雜地質(zhì)體剖分為n個(gè)四面體單元,且每個(gè)四面體密度為常數(shù),則第k個(gè)四面體在計(jì)算點(diǎn)(x,y,z)產(chǎn)生的重力異常為

(2)

其中,ρk為第k個(gè)四面體的密度值,φi為第i個(gè)三角形面的外法向量與z軸的夾角(0≤φ≤π),Si為第i個(gè)三角形的面積,Ii為第i個(gè)三角形的積分值.

(3)

為計(jì)算Ii,進(jìn)行坐標(biāo)系旋轉(zhuǎn),得到新的坐標(biāo)系(X,Y,Z),使得三角面平行于新坐標(biāo)系的XOY面,即三角面外法向量與OZ軸平行.將三角面再一次旋轉(zhuǎn),使得三角形第j邊外法向量與Y方向一致,旋轉(zhuǎn)角度φ為該邊外法向量在XOY面投影與Y軸的夾角(0≤φ≤2π).此時(shí)第i個(gè)面的積分值為

(4)

其中,Jj為第k個(gè)四面體第i個(gè)三角形第j邊的積分值.

(5)

則單個(gè)四面體重力異常Δgk的表達(dá)式為

(6)

則以非結(jié)構(gòu)網(wǎng)格四面體剖分的復(fù)雜地質(zhì)體的重力異常Δg的表達(dá)式為

(7)

2 反演方法

2.1 模型剖分

本文采用基于Delaunay三角剖分的非結(jié)構(gòu)網(wǎng)格剖分(Lawson,1986).Delaunay三角剖分依據(jù)的準(zhǔn)則是空外接圓(球)準(zhǔn)則,即在單個(gè)四面體的外接球內(nèi)不包含除本身的第5個(gè)點(diǎn).近年來(lái),三維Delaunay四面體網(wǎng)格自動(dòng)生成方法的研究也取得了一定程度的進(jìn)展,對(duì)網(wǎng)格質(zhì)量的改進(jìn),提高算法的計(jì)算效率等方面的研究也越來(lái)越多(Joe,1995; 曹書錦等,2010).為了能簡(jiǎn)單、快速地將復(fù)雜的模型剖分為規(guī)則的四面體,本文采用開源軟件TetGen,它是一款優(yōu)良的四面體網(wǎng)格生成軟件,可以對(duì)復(fù)雜的形體或地形進(jìn)行剖分(Si,2015).Tetgen剖分的輸出文件中每個(gè)四面體的三角面的頂點(diǎn)坐標(biāo)在數(shù)組中的順序不滿足正演理論的排序要求,會(huì)影響到線積分的結(jié)果.因此,需要對(duì)頂點(diǎn)進(jìn)行排序,使得從外表面方向觀察,頂點(diǎn)坐標(biāo)在數(shù)組中必須按照順時(shí)針或者逆時(shí)針排序,且各個(gè)面的排序方向相同.

目前反演分辨率要求越來(lái)越高,剖分區(qū)域的密集程度逐漸增加,使得網(wǎng)格剖分單元數(shù)迅速增加,內(nèi)存占用急劇增加,計(jì)算效率也迅速下降.為了克服這一問(wèn)題,可以通過(guò)垂向?qū)?shù)(VDR)的邊緣識(shí)別方法.垂向?qū)?shù)方法是邊界識(shí)別中應(yīng)用較為廣泛的一種方法,可以在平面上將可能存在的異常地質(zhì)體區(qū)域圈定出來(lái).將圈定的異常源區(qū)域作為靶區(qū)進(jìn)行反演,認(rèn)為其他區(qū)域不存在異常源從而不參與反演.經(jīng)過(guò)模型降維處理,可以有效地解決運(yùn)算時(shí)間長(zhǎng),內(nèi)存不足的問(wèn)題.

相關(guān)系數(shù)也是模型降維的有效方法,其表明了地下某個(gè)質(zhì)點(diǎn)引起異常的概率大小,是基于相關(guān)成像方法得到的.由重力異常的公式(1),借鑒Abdelrahman等(1989)的相鄰最小二乘剩余異常的歸一化互相關(guān)公式,則測(cè)區(qū)上實(shí)測(cè)重力異常與第q個(gè)點(diǎn)質(zhì)量在測(cè)區(qū)上的重力異常的歸一化互相關(guān)公式(郭良輝等,2009):

(8)

其中,(xi,yi,zi)為第i個(gè)觀測(cè)點(diǎn)的坐標(biāo),Δg(xi,yi,zi)為該觀測(cè)點(diǎn)的實(shí)測(cè)重力異常,Δgq(xi,yi,zi)為地下第q個(gè)點(diǎn)質(zhì)量在此觀測(cè)點(diǎn)的重力異常,N為觀測(cè)點(diǎn)總數(shù).

設(shè)剩余密度ρ>0,將(1)式代入(8)式可得

(9)

對(duì)地下空間進(jìn)行網(wǎng)格剖分,此時(shí)相關(guān)系數(shù)Cq代表了實(shí)測(cè)重力異常與地下網(wǎng)格節(jié)點(diǎn)的質(zhì)量基函數(shù)的相關(guān)性,取值范圍[-1,1],絕對(duì)值越大,表明網(wǎng)格節(jié)點(diǎn)處存在異常的概率越大.在計(jì)算所以網(wǎng)格節(jié)點(diǎn)的相關(guān)系數(shù)后,剔除絕對(duì)值相對(duì)較低的網(wǎng)格單元,即可達(dá)到模型降維的目的.

2.2 建立目標(biāo)函數(shù)

基于Tikhonov正則化反演方法,構(gòu)建目標(biāo)函數(shù):

φ=φd+λφm.

(10)

φd=‖Wd(Gm-d)‖2為數(shù)據(jù)約束項(xiàng),Wd為數(shù)據(jù)加權(quán)對(duì)角矩陣;φm=‖Wmm‖2為模型約束項(xiàng),Wm為模型加權(quán)對(duì)角矩陣;λ為正則化參數(shù),控制著數(shù)據(jù)擬合以及模型約束項(xiàng)在反演計(jì)算過(guò)程中的相對(duì)重要程度.一般使用多次試驗(yàn)、廣義交叉驗(yàn)證法(Golub et al.,1979; Li and Oldenburg,2003)或L曲線法(Hansen,1992)對(duì)正則化參數(shù)進(jìn)行取值.本文正則化參數(shù)采用多次試驗(yàn),根據(jù)反演效果確定.

在正則化反演過(guò)程中,對(duì)約束的選擇和正則化參數(shù)的選取會(huì)直接地體現(xiàn)反演結(jié)果的好壞.從目標(biāo)函數(shù)上看,數(shù)據(jù)約束項(xiàng)與模型約束項(xiàng)兩者的物理量綱不一致,正則化參數(shù)λ很難確定,本文采用最大幅值差進(jìn)行歸一化處理.

γd=dmax-dmin,

(11)

γm=mmax-mmin.

(12)

γd為數(shù)據(jù)約束項(xiàng)最大幅值差歸一化因子,d為觀測(cè)數(shù)據(jù),γm為模型約束項(xiàng)最大幅值差歸一化因子,m為密度值.此時(shí)目標(biāo)函數(shù)為

φ=‖Wd(Gm-d)/γd‖2+λ‖Wmm/γm‖2.

(13)

重力數(shù)據(jù)反演中,核函數(shù)G會(huì)隨垂向深度的增大迅速衰減,隨后逐漸平緩直至趨于零.這種現(xiàn)象在反演過(guò)程中會(huì)導(dǎo)致對(duì)模型空間淺部模型參數(shù)的權(quán)重明顯高于深部模型參數(shù)的權(quán)重,直接利用目標(biāo)函數(shù)公式進(jìn)行求解會(huì)導(dǎo)致大部分物性異常集中于地表附近,不能夠準(zhǔn)確地反映地下空間真實(shí)的物性分布狀態(tài).為了克服上述問(wèn)題,前人研究提出對(duì)反演模型參數(shù)引入深度加權(quán)函數(shù)(Li and Oldenburg,1998,2003).深度約束能有效克服重磁場(chǎng)隨深度呈指數(shù)衰減的特征.深度加權(quán)函數(shù)為

(14)

z為剖分單元的中心點(diǎn)深度,本文中取四面體各頂點(diǎn)深度的平均值;z0和β為常數(shù),z0取決于觀測(cè)面的高度,β為深度加權(quán)因子,重力反演中核函數(shù)衰減指數(shù)約為1.0～2.0.將深度加權(quán)約束加入目標(biāo)函數(shù)后形式為

φ=‖Wd(Gm-d)/γd‖2+λ‖WmWzm/γm‖2.(15)

對(duì)應(yīng)的求解方程為

(16)

在求解方程的過(guò)程中,重力反演得到的物性結(jié)果很可能不具有實(shí)際物理意義,給后續(xù)解釋帶來(lái)困難,可以在反演過(guò)程中引入物性范圍對(duì)反演解進(jìn)行物理意義上的約束.主要的物性范圍約束方式包括轉(zhuǎn)換函數(shù)法(Commer,2011)、對(duì)數(shù)轉(zhuǎn)換法(Li and Oldenburg,1996)、罰函數(shù)法(Li and Oldenburg,2003).本文使用強(qiáng)制約束的方法,將反演密度值約束在給定的范圍,其上下限的范圍根據(jù)實(shí)際地質(zhì)資料的密度范圍確定,在求解的每一次迭代過(guò)程中,判斷求解的密度值是否超過(guò)了理論的界限,若超過(guò)則賦值為界限值,再將修改過(guò)的解方程作為下一次的迭代初值,直到符合要求.

共軛梯度法具有計(jì)算速度快、內(nèi)存需求少,能夠避免大型矩陣的乘積計(jì)算的特點(diǎn),是求解線性方程組的有效方法之一(Wilson et al.,2011; 陳少華等,2013; Hou et al.,2015).為了改善迭代法的收斂性,增強(qiáng)解的穩(wěn)定性,求解方程組的方法選用預(yù)條件共軛梯度法(PCG),將原方程Ax=b替換為SAx=Sb,用深度加權(quán)函數(shù)即式w(z)作為預(yù)條件因子S.

另外,為了進(jìn)一步提高反演的效率,我們采用并行進(jìn)行計(jì)算加速.OpenMP是共享存儲(chǔ)的并行編程模型,具有可移植和易編程性等特點(diǎn),適用于循環(huán)級(jí)并行等任務(wù),對(duì)程序代碼改動(dòng)小,在科學(xué)計(jì)算上占據(jù)優(yōu)勢(shì).因此本文采用OpenMP并行計(jì)算方法來(lái)提升三維密度的反演效率.

3 模型試驗(yàn)

3.1 單個(gè)模型

對(duì)于一個(gè)邊長(zhǎng)200 m×200 m×200 m的長(zhǎng)方體,參數(shù)如Li和Oldenburg(1998),平面測(cè)網(wǎng)在x、y方向?yàn)?～1000 m,點(diǎn)距為25 m,高度位于-100 m,將模型頂面深度置于50 m,x、y方向置于400～600 m處,模型(圖2a)的剩余密度為1.0 g·cm-3.地區(qū)的深度為500 m,對(duì)整個(gè)地區(qū)進(jìn)行均勻剖分,剖分單元為8926個(gè),其剖面示意圖為圖2b,除模型外其他區(qū)域的剩余密度為0.0 g·cm-3,計(jì)算重力異常(圖2c).為了驗(yàn)證四面體公式的準(zhǔn)確性,在相同參數(shù)下利用直立長(zhǎng)方體的解析式計(jì)算重力異常(圖2d),兩者的殘差結(jié)果表明(圖2e),利用四面體公式計(jì)算的重力異常與解析式的計(jì)算結(jié)果一致.

圖2 直立長(zhǎng)方體模型正演

根據(jù)正演的重力異常數(shù)據(jù)作為觀測(cè)數(shù)據(jù)進(jìn)行反演模型試驗(yàn),首先不添加任何約束,僅使用最小長(zhǎng)度作為模型的約束項(xiàng),利用預(yù)條件共軛梯度法求解方程,正則化參數(shù)λ為1.0,反演結(jié)果的垂直切片如圖3a.為了更直觀地分析,在x=500 m處對(duì)模型切片(圖3b),異常偏離了原本位置,趨近于地表,這是由于正演算子核函數(shù)G會(huì)隨著垂向深度的增大而衰減,因此需要添加深度約束校正.

圖3 最小長(zhǎng)度反演結(jié)果

為了確定正則化參數(shù),在添加深度約束前進(jìn)行歸一化,使兩者保持在同一量級(jí).深度約束根據(jù)(14)式,設(shè)置參數(shù)z0為15,β為2.0,正則化參數(shù)λ為1.0,反演結(jié)果的垂直切片如圖4a.分別取深度為50 m,150 m,250 m,350 m的水平切片組合(圖4b)以及沿x=500 m的垂直切片(圖4c),模型異常在深度上集中于150 m處,基本符合直立長(zhǎng)方體模型的大致范圍.但反演結(jié)果存在拖尾現(xiàn)象,模型下方反演密度值逐漸擴(kuò)散并減小.因?yàn)榈叵驴臻g所有剖分單元在模型反演過(guò)程中都提供了貢獻(xiàn),因此預(yù)測(cè)的密度值會(huì)依據(jù)貢獻(xiàn)大小分散于整個(gè)地下空間,模型的密度也會(huì)明顯偏小,Li和Oldenburg(1998)的模型試驗(yàn)的反演結(jié)果同樣低于理論值.

圖4 添加深度加權(quán)約束反演結(jié)果

在添加了深度加權(quán)約束后,模型處的異常集中于0.5～0.7 g·cm-3處,但地下區(qū)域的邊緣位置大多為負(fù)值,也有個(gè)別較大的極值點(diǎn),為了限制物性的范圍,添加物性約束,將地下區(qū)域反演結(jié)果強(qiáng)制限定在0.0～1.0 g·cm-3之間,反演三維切片如圖5a.在與圖4相同位置上進(jìn)行水平切片和沿x=500 m的垂直切片(圖5b、圖5c),異常依然集中于直立長(zhǎng)方體范圍內(nèi),但由于進(jìn)行了物性約束,邊緣位置的密度值接近0.0 g·cm-3,在反演中的貢獻(xiàn)減小,所以模型周圍的密度反演值相對(duì)于圖4有所收斂,解決了在深度方向上的拖尾現(xiàn)象,模型中心的異常值也相對(duì)更加接近理論值.

圖5 添加物性約束后反演結(jié)果

正則化反演的過(guò)程中,構(gòu)造核函數(shù)、求解方程組占用了大部分的時(shí)間,隨著地下剖分單元的增加,反演越來(lái)越長(zhǎng),且會(huì)產(chǎn)生內(nèi)存不足的問(wèn)題,解決這類問(wèn)題可以使用模型降維的方法.首先通過(guò)垂向?qū)?shù),在平面內(nèi)大致圈定出異常存在的區(qū)域,選擇大于某一常數(shù)值的區(qū)域作為模型可能存在的區(qū)域,區(qū)域外的四面體單元不參與反演,利用剩下的單元進(jìn)行反演.反演的異常體(圖6)相較于圖5更加明顯,同時(shí)矩陣內(nèi)單元數(shù)約為原來(lái)的四分之一,大大節(jié)省了內(nèi)存空間.根據(jù)反演的密度預(yù)測(cè)值進(jìn)行重力異常的數(shù)據(jù)重建(圖6d),與觀測(cè)得重力異常(圖2c)相比,殘差(圖6e)最大值為 0.096 mGal, 最小值為0.012 mGal, 平均值為0.048 mGal, 均方根為0.054 mGal.重力異常的重建效果較好.下文中模型試驗(yàn)反演均對(duì)重力數(shù)據(jù)進(jìn)行了重建,效果較好,但限于篇幅不再展示.

圖6 垂向一階導(dǎo)數(shù)降維反演結(jié)果

為了測(cè)試深度加權(quán)函數(shù)對(duì)不同深度模型的反演效果,將同樣的直立長(zhǎng)方體模型保持平面位置不變,頂面埋深分別改為100 m和150 m.對(duì)于埋深為100 m的模型,設(shè)置參數(shù)z0為15,β為1.5,正則化參數(shù)λ為1.0,反演三維結(jié)果為圖7a,根據(jù)水平和x方向的切片顯示(圖7b、圖7c),異常集中于模型位置處; 對(duì)于埋深為150 m的模型,設(shè)置參數(shù)z0為5,β為2.0,λ為10.0,反演三維結(jié)果為圖7d—f,異常中心位于模型位置處.從頂面埋深50 m到100 m,再到150 m,可以看到需要的深度加權(quán)越來(lái)越強(qiáng),但反演的效果隨著深度變化越來(lái)越弱,異常中心逐漸向上偏移,在埋深150 m處除了要加強(qiáng)深度加權(quán)函數(shù),還需要放大正則化參數(shù),加強(qiáng)模型項(xiàng)的影響,但效果依然較為發(fā)散.

圖7 不同埋深的直立長(zhǎng)方體反演

除了垂向?qū)?shù),還可以利用正演數(shù)據(jù)和剖分單元的重心坐標(biāo)計(jì)算地下區(qū)域每個(gè)單元的相關(guān)系數(shù),根據(jù)相關(guān)系數(shù)選取四面體單元,將相關(guān)系數(shù)低的單元從矩陣中排除達(dá)到降維的效果.以上述所展示的三個(gè)不同埋深的直立長(zhǎng)方體為例,利用相關(guān)系數(shù)進(jìn)行降維,不改變參數(shù)值,反演結(jié)果如圖8所示.埋深50 m的模型反演結(jié)果(圖8a—c)顯示,在深度上密度值集中于模型范圍內(nèi),但偏離模型中心,向下偏移,下方有拖尾現(xiàn)象.埋深100 m的反演結(jié)果(圖8d—f)和埋深150 m的模型反演結(jié)果(圖8g—i),異常基本集中于模型的中心,結(jié)果收斂.相較于利用垂向?qū)?shù)降維后的深度約束反演結(jié)果,引入相關(guān)系數(shù)降維,反演效果在埋深較深的情況下取得了更理想的效果,形態(tài)更加收斂,但對(duì)于埋深較淺的模型,效果向下偏移,形態(tài)發(fā)散.因此,針對(duì)埋深較淺的異常體,利用垂向?qū)?shù)降維的深度加權(quán)反演結(jié)果更好; 針對(duì)埋深較深的異常體,可以通過(guò)相關(guān)系數(shù)降維的深度加權(quán)得到理想反演效果.

圖8 不同埋深的相關(guān)系數(shù)降維反演結(jié)果

3.2 組合模型

實(shí)際地下區(qū)域的異常會(huì)更加復(fù)雜,為了更加符合現(xiàn)實(shí)情況,構(gòu)建多個(gè)模型組成復(fù)雜的異常體.在地下中心位置(500,400,200)處建立一個(gè)直立長(zhǎng)方體,邊長(zhǎng)均為200 m,以中心位置(500,700,150)為球心建立一個(gè)球體,半徑為100 m,此組合模型的正演數(shù)據(jù)(圖9a)為觀測(cè)數(shù)據(jù)進(jìn)行反演,分別采用垂向一階導(dǎo)數(shù)(圖9c—e)和相關(guān)系數(shù)降維(圖9f—h)兩種方法.反演結(jié)果顯示,由于兩個(gè)模型的埋深較淺,在縱向上使用垂向?qū)?shù)的效果更好,但在橫向上相對(duì)發(fā)散,區(qū)分效果相對(duì)一般; 在縱向上用相關(guān)系數(shù)降維的效果基本符合實(shí)際位置,但存在拖尾的現(xiàn)象,而異常在橫向收斂于模型范圍內(nèi),可以更加清楚地分辨兩個(gè)模型的位置.為了測(cè)試噪聲對(duì)反演的影響,向觀測(cè)數(shù)據(jù)中加入信噪比為5%的高斯白噪聲(圖9b),分別采用垂向一階導(dǎo)數(shù)(圖9i—k)和相關(guān)系數(shù)降維(圖9l—n)進(jìn)行三維反演,密度相對(duì)于不加噪的結(jié)果略微分散,異常中心最大值較低,但依舊符合實(shí)際位置,因此本文所提方法具有較好的穩(wěn)定性.

圖9 組合模型不同降維方法的反演結(jié)果

3.3 起伏地形下模型

地形的起伏對(duì)三維密度反演效果存在一定影響,將地形因素考慮到模型建模中.使用函數(shù)模擬小幅度隆起的地形(圖10a),根據(jù)地形構(gòu)建模型(圖10b),剖分為9906個(gè)單元,觀測(cè)面、模型等數(shù)據(jù)與前文保持一致.

圖10 起伏地形模擬和模型構(gòu)建

模型選用3.2節(jié)的組合模型,垂向?qū)?shù)降維的結(jié)果異常集中于模型的范圍內(nèi)(圖11a—c),在深度上結(jié)果較為收斂,但由于兩個(gè)模型的位置較為接近,反演效果在橫向上區(qū)分度一般,模型中間的區(qū)域與異常大小相近; 相關(guān)系數(shù)降維的反演結(jié)果同樣基本符合模型位置(圖11d—f),在深度上存在拖尾,但在橫向上,能夠明顯地區(qū)分兩個(gè)異常的中心位置.因此,本文所使用的方法同樣適用于起伏地形.

圖11 起伏地形下組合模型的反演結(jié)果

4 實(shí)際資料反演

氧化鐵礦床是非常理想的地球物理目標(biāo),尤其是在區(qū)域尺度上的重磁測(cè)量中(Ives and Mickus,2019; Vatankhah et al.,2022).本文選取美國(guó)密蘇里州東南部以Pea Ridge為中心的氧化鐵礦床,面積為35×37 km2,包括Kratz Spring和Bourbon礦床.2014年美國(guó)地質(zhì)調(diào)查局USGS(McCafferty,2016)聯(lián)合CGG利用HeliFALCON AGG系統(tǒng)進(jìn)行的航空重力梯度測(cè)量中,測(cè)區(qū)內(nèi)飛行總長(zhǎng)度為3537.7 km,測(cè)線線距為400 m,地形數(shù)據(jù)(圖12a)由兩臺(tái)激光掃描儀繪制,航空重力測(cè)量高度(圖12b)由雷達(dá)高度計(jì)以0.1 s為間隔測(cè)量,距離地面90 m,通過(guò)測(cè)得的兩個(gè)重力梯度的水平曲率分量計(jì)算測(cè)區(qū)的重力梯度.根據(jù)重力梯度數(shù)據(jù)利用傅里葉域變換方法計(jì)算地形校正密度為2.67 g·cm-3的布格重力異常數(shù)據(jù)(圖12c).

圖12 密蘇里州東南部礦區(qū)地形及重力異常數(shù)據(jù)

Kratz Spring和Bourbon地區(qū)的隱伏礦床全部被古生代地層覆蓋(Day et al.,2016),根據(jù)物性測(cè)量結(jié)果得到基底起始密度: 2.78 g·cm-3(Kratz Spring)和2.83 g·cm-3(Bourbon).Kratz Spring為氧化鐵-磷灰石(IOA)型礦床,礦床的鐵氧化物體積較小,與周邊礦區(qū)鐵礦相比,特點(diǎn)是密度低,磁性巖石較少.在礦床南部有一個(gè)大而連續(xù)的、向南傾斜的致密源區(qū)域,可能為一個(gè)以赤鐵礦為主的潛伏礦,該巖性密度最大的部分大約在古生代接觸面以下300 m處,延伸到海平面以下超過(guò)2000 m的深度.Bourbon為氧化鐵-銅-金(IOCG)類型的礦床,形成于地表以下400 m至490 m處的前寒武紀(jì)/古生代界面,是密蘇里州東南部規(guī)模最大的礦床.礦床頂部以磁鐵礦和赤鐵礦為主,呈碗狀的大體積高密度高磁化率巖石,直徑為4200 m,深至海平面以下1200 m(McCafferty et al.,2016).

將反演的主要目標(biāo)定為Kratz Spring和Bourbon礦床,根據(jù)前人資料顯示,兩者密度變化較為顯著的深度為500～1200 m左右,深度密度變化不大,因此將反演地區(qū)的深度設(shè)置為2000 m,根據(jù)地形數(shù)據(jù)組合為反演區(qū)域,并對(duì)該區(qū)域進(jìn)行四面體剖分,剖分個(gè)數(shù)為39042個(gè),觀測(cè)面x方向從650～686 km,y方向從4206～4242.8 km,點(diǎn)距為400 m,總點(diǎn)數(shù)8463個(gè).通過(guò)垂向一階導(dǎo)數(shù)選取重點(diǎn)區(qū)域,降維后對(duì)礦區(qū)進(jìn)行三維重力密度反演,深度加權(quán)的參數(shù)z0為10,參數(shù)β為 2.0,正則化參數(shù)λ為2.0,物性約束范圍為-0.3～0.3 g·cm-3.圖13和圖14分別為采用垂向?qū)?shù)降維和相關(guān)系數(shù)降維的反演結(jié)果及重力數(shù)據(jù)重建.兩種降維方法重建的重力異常基本與觀測(cè)數(shù)據(jù)特征一致,僅在實(shí)測(cè)范圍的邊緣部分存在一些差異,重力數(shù)據(jù)重建的均方差分別為0.803 mGal和0.351 mGal(圖13a、圖13b; 圖14a、圖14b),顯示出了良好的重力數(shù)據(jù)擬合效果.

圖13 密蘇里州東南部三維反演(垂向一階導(dǎo)數(shù)降維)

圖14 密蘇里州東南部三維反演(相關(guān)系數(shù)降維)

兩種不同降維方法的三維反演結(jié)果均顯示了密蘇里東南部存在的高密度巖體(圖13c; 圖14c),其形態(tài)基本一致,在500 m、1000 m、1500 m處的水平切片(圖13d—f; 圖14d—f)也顯示出同樣的高密度巖體分布特征.兩種不同降維方法反演結(jié)果的垂直切片顯示(位置如圖12c): 在過(guò)Kratz Spring的剖面兩者的反演結(jié)果較為一致(圖15a、圖15c),并與前人的反演結(jié)果(McCafferty et al.,2016)比較形態(tài)基本吻合.其北側(cè)沉積層下存在較為明顯的中密度體,南側(cè)的淺部為中密度體,延伸至深部為高密度的巖體,且有向南傾斜的趨勢(shì).而在過(guò)Bourbon的剖面兩者的反演結(jié)果有所不同(圖15b、圖15d),其中采用垂向?qū)?shù)降維反演的高密度巖體集中在0～1000 m深度范圍內(nèi),形狀呈類似的碗狀,與前人的反演結(jié)果(McCafferty et al.,2016)相比較為吻合,形狀略有不同.而通過(guò)相關(guān)系數(shù)降維的反演受相關(guān)系數(shù)影響,礦體的反演結(jié)果偏深部.

圖15 Kratz spring 和 Bourbon 反演結(jié)果

兩種方法的反演效果略有區(qū)別,主要原因是由于模型降維的篩選方式不同: 垂向一階導(dǎo)數(shù)常用于邊界識(shí)別,平面上可以清晰反應(yīng)異常體位置,篩選方式是根據(jù)垂向一階導(dǎo)數(shù)值保留垂直方向上的所有四面體,因此在邊緣位置的密度區(qū)分明顯; 而相關(guān)系數(shù)反應(yīng)了剖分單元引起異常的概率大小,計(jì)算結(jié)果在深度上集中于異常體底部,因此底面模型單元的保留個(gè)數(shù)會(huì)大于表面單元數(shù),反演結(jié)果在深部較為明顯,在淺部反演會(huì)受到影響.

5 結(jié)論

針對(duì)起伏地形和復(fù)雜地質(zhì)體的重力三維密度反演問(wèn)題,本文提出了基于Delaunay四面體非結(jié)構(gòu)網(wǎng)格剖分的重力三維密度正則化反演方法,該方法引入了模型降維、模型約束以及OpenMP并行算法等,并采用預(yù)條件共軛梯度法求解.通過(guò)模型試驗(yàn)以及對(duì)美國(guó)密蘇里東南部礦區(qū)實(shí)測(cè)資料進(jìn)行反演驗(yàn)證,得到了如下結(jié)論:

(1) 四面體非結(jié)構(gòu)剖分可以較為準(zhǔn)確地?cái)M合起伏地形,減小了地形對(duì)反演的影響.最小長(zhǎng)度約束、深度加權(quán)約束和物性范圍約束可以有效改善四面體非結(jié)構(gòu)網(wǎng)格模型的反演效果,準(zhǔn)確地反演起伏地形下的地質(zhì)體.通過(guò)深度加權(quán)約束,解決了反演結(jié)果趨于地表的現(xiàn)象,基本符合異常體的大致范圍; 通過(guò)物性范圍約束,解決了深部的拖尾現(xiàn)象,反演收斂于異常體中心.

(2) 引入垂向?qū)?shù)和相關(guān)系數(shù)的模型降維有效減小解的空間,能夠突出目標(biāo)區(qū)域,提高反演效率.對(duì)于埋深較淺的異常體,引入垂向?qū)?shù)的模型降維方法,可以更加明顯地反映地下異常體的位置,能有效地解決深度加權(quán)約束中的拖尾現(xiàn)象; 對(duì)于埋深較深的異常體,引入相關(guān)系數(shù)的模型降維方法效果更好,在橫向上對(duì)于異常體之間的區(qū)分度更好.

致謝感謝論文匿名評(píng)審專家和編輯部提出的寶貴修改意見,感謝TetGen官網(wǎng)提供的軟件,感謝美國(guó)地質(zhì)調(diào)查局(USGS)提供的重力異常數(shù)據(jù).