基于多重水平集的三維磁性體邊界反演

肖 曉 段雅婷 胡雙貴 湯井田 謝 勇 劉長生

(①中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室，湖南長沙 410083； ②有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室，湖南長沙 410083； ③中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長沙 410083； ④自然資源部覆蓋區(qū)深部資源勘查工程技術(shù)創(chuàng)新中心，湖南長沙 410083； ⑤長沙航空職業(yè)技術(shù)學(xué)院，湖南長沙 410124)

0 引言

衛(wèi)星、航空、地面、井中等多維空間的磁場數(shù)據(jù)測量技術(shù)迅速發(fā)展并逐漸成熟，配套的數(shù)據(jù)處理和解釋技術(shù)也在不斷更新和完善，這勢必拓寬磁法勘探的應(yīng)用領(lǐng)域，提高磁法的探測能力[1-4]。磁場數(shù)據(jù)三維反演不僅能夠提供磁異常體的空間位置，還能獲得磁性異常體的形狀、物性分布等定量信息，在磁場數(shù)據(jù)處理和解釋中具有重要意義，一直是該領(lǐng)域的研究熱點問題。

目前的磁異常反演包括磁化率和磁性界面兩方面的內(nèi)容。磁化率反演首先固定反演網(wǎng)格，假設(shè)反演單元內(nèi)磁化率為常數(shù)或者線性分布，建立反演目標函數(shù)，然后最小化反演目標函數(shù)，從而獲得目標的磁化率分布[5-6]。Li等[6]利用一個或多個適當?shù)募訖?quán)函數(shù)將先驗信息引入目標函數(shù)；Pilkington[7]對目標函數(shù)進行平滑及深度加權(quán);Miguel[8]通過地質(zhì)統(tǒng)計學(xué)模型實現(xiàn)了巖性約束下的重磁數(shù)據(jù)聯(lián)合反演；Fregoso等[9]實現(xiàn)了重磁數(shù)據(jù)的交叉梯度聯(lián)合反演。上述方法可定量描述地質(zhì)構(gòu)造的磁化率分布。磁性界面反演方法一般假定已知地下異常體的磁化率，通過反演得到異常體的邊界和位置。管志寧等[10]提出了頻率域磁性單界面反演方法，推導(dǎo)了常磁性與變磁性單界面反問題解的近似表達式；Pilkington等[11]提出了利用重磁資料確定地殼界面起伏形態(tài)；Wang等[12]對多面體的頂點位置進行反演；劉沈衡等[13]利用歐拉算法反演磁性界面；Fullagar等[14]通過混合參數(shù)化反演得到目標體磁化率的分布及垂直界面的位置。趙文舉等[15]通過BP神經(jīng)網(wǎng)絡(luò)預(yù)測磁性體頂面埋深；鄭強等[16]基于磁梯度張量特征值成功地識別了磁性體邊界。然而，由于缺乏地下異常體拓撲結(jié)構(gòu)和數(shù)量信息，現(xiàn)有的磁邊界反演算法無法靈活描述復(fù)雜地下磁性目標體的幾何形狀。

因此，需要一種無需過多關(guān)于異常體拓撲結(jié)構(gòu)和數(shù)量的先驗信息就可以直接反演地下磁異常體空間位置和邊界的算法?；谒郊椒╗17-21]的磁邊界反演算法正好符合這一要求，該方法僅需已知目標體的平均磁化率就可以反演得到地下異常體的邊界及幾何形狀。Isakov等[22]將水平集方法用于重力數(shù)據(jù)反演；Zheglova等[23]將水平集用于地震旅行時層析成像，實現(xiàn)了邊界的二維重建；Li等[24]在磁性數(shù)據(jù)的反演中引入水平集，用兩個水平集對地下磁性體的邊界進行了三維反演；在礦產(chǎn)勘查中，Li等[25]在地表數(shù)據(jù)與鉆孔磁資料的聯(lián)合反演中引入水平集，實現(xiàn)了目標圈定。基于水平集的反演可以得到異常體明確的邊界，將水平集引入磁邊界反演，是改善磁邊界反演效果的重要途徑。然而，目前已發(fā)表的基于水平集的磁性目標體邊界反演算法僅使用了兩個水平集，而實際地球物理勘探的目標往往是多個具有不同磁化率的磁性體。

針對上述問題，本文在Li等[24]的基礎(chǔ)上提出了一種基于多重水平集的磁邊界反演算法，該算法具有抗噪性強、精度高、多目標識別的優(yōu)點。首先基于多個水平集函數(shù)建立磁界面反演目標函數(shù)，并詳細推導(dǎo)了在多重水平集反演構(gòu)架下磁邊界反演的實現(xiàn)過程；然后采用任意四面體單元磁法解析解算法[24]進行高精度正演計算，構(gòu)建正演網(wǎng)格與反演網(wǎng)格的映射關(guān)系，并基于加權(quán)基本無振蕩(WENO)格式對水平集函數(shù)進行更新及重新初始化；最后，通過多個理論模型算例驗證本文反演算法的有效性，并分析初始模型對反演結(jié)果的影響。

1 方法與原理

1.1 基于磁化率的多重水平集表示

假設(shè)非磁性背景下，區(qū)域Ω中包含N個磁性體，且每個磁性體具有恒定的磁化率，則磁化率分布可表示為

(1)

上述模型的磁化率可用水平集公式表示為

(2)

式中：φ(·)表示水平集函數(shù)；H(·)表示Heaviside函數(shù)[25]。其表達式分別為

(3)

(4)

(5)

(6)

式中δ(·)是Dirac delta函數(shù)。

1.2 基于多重水平集的磁邊界反演算法

反演的總目標函數(shù)Ei可表示為數(shù)據(jù)擬合差泛函Ed和正則化項Er的線性組合

Et=Ed+αEr

(7)

1.2.1 數(shù)據(jù)失配項及其導(dǎo)數(shù)

(8)

式中：M是觀測點數(shù)；σk是與第k個數(shù)據(jù)相關(guān)的誤差標準偏差，數(shù)據(jù)擬合差的最優(yōu)值一般為0.5。

(9)

(10)

(11)

(12)

(13)

令

(14)

則

(15)

1.2.2 正則化項及其導(dǎo)數(shù)

為了進一步降低反演的多解性，本文在水平集反演目標函數(shù)中引入Tikhonov正則化項

(16)

式中P代表水平集個數(shù)。Er的Fréchet導(dǎo)數(shù)為

(17)

(18)

1.3 水平集函數(shù)的更新及重新初始化

對均勻網(wǎng)格，利用WENO格式對水平集方程和重新初始化方程進行空間離散[29-34]，并利用三階精確總變差遞減龍格—庫塔格式(TVD-RK3)[32]進行時間離散化。

1.3.1 水平集方程的近似

對于式(18)，使用TVD-RK3格式及Godunov法對H(φi)=Fi|φi|進行空間離散化，得到

(19)

式中：

使用TVD-RK3格式進行時間離散化。例如,對水平集方程φt+Fi|φ|=0，首先對φ(l)執(zhí)行一次歐拉計算，得到t(l+1)時刻的解

(20)

(21)

(22)

通過線性平均更新φ(l+1)

(23)

并建立一個Courant-Friedrichs-Lewy(CFL)條件[34]

(24)

保持更新的穩(wěn)定性。在本文反演中，正則化參數(shù)α遠小于max|Fi|，所以設(shè)

(25)

1.3.2 重新初始化方程的近似

一般情況下，水平集函數(shù)在通過式(18)進行迭代時，并不保留其符號距離屬性，所以Sussman等[19]引入了重新初始化方程

(26)

(27)

2 反演算例

2.1 反演流程

本文采用任意四面體單元磁法解析解算法[24]進行高精度正演計算，利用正、反演網(wǎng)格之間的物性映射實現(xiàn)正演網(wǎng)格與反演網(wǎng)格相互獨立運行，并基于WENO格式對水平集函數(shù)進行更新及重新初始化，具體反演流程如下。

(1)初始化水平集函數(shù)φi；

(4)基于WENO格式計算總目標函數(shù)Et的負梯度方向-?Et/?φi，并更新水平集函數(shù)φi；

(5)重新初始化水平集函數(shù)φi，使其保持帶符號距離屬性；

(6)返回步驟(2)，直至迭代收斂或達到最大迭代次數(shù)。

2.2 兩個水平集模型

設(shè)計一個由具有不同磁化率的兩個傾斜棱鏡體組成的模型，分析基于兩個水平集公式的磁邊界反演效果。首先根據(jù)先驗信息確定兩個不同的磁化率，并使用基于兩個水平集公式的磁邊界反演算法確定異常體數(shù)目并恢復(fù)磁性體邊界。

設(shè)計模型如圖1a所示，由κ1=0.04SI和κ2=0.08SI的兩個傾斜棱鏡體組成。圖1b為加噪總磁異常圖。初始猜測如圖1c所示。反演中引入兩個水平集函數(shù)φ1和φ2，初始猜測分別為

(28)

φ2=1-

(29)

圖1d是使用多水平集反演得到的模型，其數(shù)據(jù)擬合差Ed=0.5022，可見反演結(jié)果與原模型的形狀和位置吻合較好。圖1e是磁場強度B0=50000nT、磁傾角為75°、磁偏角為25°的背景感應(yīng)場下，反演模型(圖1d)在觀測面正演的磁異常圖。圖1f為反演模型的正演磁異常(圖1e)與原模型加噪磁異常(圖1b)之差。反演結(jié)果的更多細節(jié)見圖2。顯然，使用兩個水平集的反演結(jié)果很好地恢復(fù)了每個棱鏡體的傾角、空間位置及形態(tài)，且反演誤差與預(yù)期的一樣，呈隨機分布特征。

圖1 兩個不同磁化率的傾斜棱鏡模型及兩個水平集計算結(jié)果(a)模型示意圖； (b)加噪總磁異常圖； (c)初始猜測模型； (d)兩個水平集的反演模型； (e)圖d正演磁異常圖； (f)圖e與圖b的差值

圖2 兩個不同磁化率的傾斜棱鏡體模型使用兩個水平集的反演結(jié)果細節(jié)展示(a)y=300m截面； (b)y=700m截面； (c)z=150m平面圖； (d)z=250m平面圖虛線表示模型的真實邊界，實線表示模型的恢復(fù)邊界

將本文反演結(jié)果與Li等[24]的反演結(jié)果做對比，見圖3?？梢钥闯觯疚牟捎没谌我馑拿骟w單元的磁法解析解算法[26]進行高精度正演，并引入HJ-WENO格式對水平集方程進行更新及重新初始化，反演得到的磁異常體邊界與模型吻合更好。

圖3 本文反演結(jié)果與Li等[24]的反演結(jié)果對比(a)y=300m截面(過磁性體①)； (b)y=700m截面(過磁性體②)虛線為模型的輪廓線，實線為本文反演磁性體邊界，*線表示Li等反演的磁性體邊界

2.3 三個水平集模型

實際地下地質(zhì)情況是非常復(fù)雜的，需分析多個磁性體存在的情況下，本文方法是否有效。對于具有不同磁化率值的磁性體，需要引入更多的水平集函數(shù)，對其進行磁邊界反演。

假設(shè)模型由三個具有不同磁化率的直立長方體組成(圖4a)，且磁化率值已知。圖4b為加噪總磁異常圖。根據(jù)圖4b，在反演算法中引入三個水平集函數(shù)φ1、φ2和φ3，初始猜測(圖4c)分別為

φ1=1-

(30)

φ2=1-

(31)

φ3=1-

(32)

圖4 三個不同磁化率模型及三個水平集計算結(jié)果(a)模型示意圖； (b)加噪總磁異常圖； (c)初始猜測模型； (d)三個水平集的反演模型；(e)圖d正演磁異常圖； (f)圖e與圖b的差值圖a中磁性體①、②、③的長方體中心坐標分別為(500m,150m,175m)、(500m,475m,350m)和(500m,825m,300m)，大小分別為400m×200m×200m、200m×150m×200m和400m×150m×200m，磁化率分別為0.04、0.08、0.12SI

圖5 圖4d細節(jié)展示(a)y=200m截面； (b)y=500m截面； (c)y=800m截面； (d)z=100m平面圖虛線為真實模型的輪廓線，實線為本文反演結(jié)果

2.4 深部異常體模型

前面兩個模型中的磁異常體深度較小，說明了本反演算法對淺層磁性礦藏具有很好的圈定作用。為了進一步驗證該方法對深部目標體邊界的圈定效果，設(shè)計了圖6a所示的模型。

模型包括兩個磁性長方體，其中心坐標分別為 (2000m,1500m,900m)、(2950m,2550m,825m)，大小分別為1000m×500m×400m、700m×300m×350m,極化率已知，分別為0.16、0.30SI。基于加噪總磁異常數(shù)據(jù)(圖6b)，反演中引入兩個水平集函數(shù)φ1和φ2，其初始猜測(圖6c)分別為

φ1=1-

(33)

φ2=1-

(34)

圖6d為反演結(jié)果，其數(shù)據(jù)擬合差Ed=0.5916。圖6e和圖6f分別為反演模型的正演磁異常圖及誤差。圖7是圖6d的細節(jié)展示。由圖可見，反演結(jié)果在各個方向上對目標體的邊界都吻合得很好，展現(xiàn)了真實磁性體的形狀和空間位置。

圖6 深部磁異常體模型及反演結(jié)果(a)模型示意圖； (b)加噪總磁異常圖； (c)初始猜測模型； (d)基于兩個水平集的反演結(jié)果； (e) 圖d正演磁異常圖； (f)圖e與圖b的差值

圖7 圖6d細節(jié)展示(a)z=800m平面圖； (b)x=2000m截面；(c)x=3000m截面圖。虛線表示原始模型輪廓線，實線表示反演模型邊界

本算例驗證了本文反演算法對深部礦產(chǎn)的圈定效果很好。

2.5 水平集數(shù)目對反演結(jié)果的影響

基于多重水平集方法的磁界面反演算法的關(guān)鍵假設(shè)是地下異常體的磁化率已知。實際反演中，一般可以基于已有的地質(zhì)知識或其他先驗信息得到研究區(qū)磁性體確切的磁化率值，但無法確切知道某一區(qū)域具體對應(yīng)哪個磁化率值。復(fù)雜的地質(zhì)環(huán)境中，這個不可預(yù)知的因素可能會極大地限制多重水平集方法的應(yīng)用。本節(jié)針對圖4a所示模型，假設(shè)僅知道一個磁化率值，分析水平集數(shù)目的選擇對反演結(jié)果的影響。

假設(shè)圖4a中僅知道②號異常體的磁化率(κ=0.08SI)，并據(jù)此形成初始猜測(圖8a)，并在反演中引入一個水平集函數(shù)

φ=1-

(35)

圖8b為引入單水平集的反演結(jié)果。可見即使初始模型只假設(shè)了一個磁源體，反演結(jié)果也能準確地反映有三個獨立的磁源體，說明水平集的個數(shù)對反演結(jié)果影響不大。

圖9是圖8b的細節(jié)展示。對比圖9a與圖6a可以清晰地看出，當κ=0.04SI的①號長方體被賦予大于其實際磁化率值(κ=0.08SI)時，反演得到的磁性邊界范圍更小、深度更大；對比圖9b與圖6a可以看出，當磁化率為κ=0.12SI的②號長方體被賦予低于其實際磁化率值(κ=0.08SI)時，反演得到的邊界范圍更大，且深度小于實際深度；圖9c中，反演的③號磁性體邊界與模型(圖6a)幾乎完全一致，這是因為反演時設(shè)定的磁化率值與實際情況完全一致。

圖8 圖6a模型單水平集反演結(jié)果(a)初始猜測； (b)反演模型

圖9 圖8b細節(jié)展示(a)y=200m截面(過①號磁性體)； (b)y=500m截面(過②號磁性體)； (c)y=800m截面(過③號磁性體)虛線表示原模型磁性體輪廓線，實線表示反演模型邊界

該模型反演結(jié)果說明，基于水平集方法反演時模型的磁化率與真實值差別太大時，恢復(fù)的磁性體邊界會出現(xiàn)較大偏差：如果反演模型的磁化率賦值低于實際值，反演結(jié)果中磁性源的體積會大于實際體積，且位置偏深；反之則會導(dǎo)致反演的磁源體積小于實際體積，且位置偏淺。

由本節(jié)反演結(jié)果可知，反演初始模型的各個子域磁化率值是否準確，對反演結(jié)果影響較大，而反演初始模型水平集的數(shù)目對反演結(jié)果的影響較小。因此，在實際應(yīng)用中，需要利用先驗信息盡量準確地估計各個區(qū)域的磁化率值，以各區(qū)域的平均磁化率作為各子域的磁化率進行多水平集反演，才可得到比較可靠的磁邊界位置。

3 結(jié)論

本文基于多重水平集的磁界面反演算法，實現(xiàn)了對磁化率已知的磁性目標體空間位置和幾何形狀的反演。

(1)本文算法中，磁化率值由先驗信息確定，采用任意四面體單元磁法解析解算法[24]進行高精度正演計算，并基于WENO格式對水平集函數(shù)進行更新及重新初始化。

(2)具有兩個水平集和三個水平集的反演算例驗證了本算法的有效性與可靠性；深部異常體模型反演算例驗證了本算法對深部礦產(chǎn)資源的圈定具有良好效果。

(3)針對實際生產(chǎn)中，難以準確獲取地下目標磁化率的情況，討論了水平集數(shù)目的確定對反演結(jié)果的影響。即便水平集數(shù)目與實際不符，反演結(jié)果仍能準確劃分磁性體區(qū)域，但各反演磁性體的深度和邊界存在一定誤差。若子域被賦予偏大的磁化率值時，反演得到的磁性體范圍小于實際情況，且深度偏小；反之，若子域被賦予較小磁化率值時，則反演磁性體的范圍大于實際情況，且深度偏大。