肖 曉 段雅婷 胡雙貴 湯井田 謝 勇 劉長(zhǎng)生

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

0 引言

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

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

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

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

1 方法與原理

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

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

(1)

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

(2)

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

(3)

(4)

(5)

(6)

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

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

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

Et=Ed+αEr

(7)

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

(8)

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

(9)

(10)

(11)

(12)

(13)

令

(14)

則

(15)

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

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

(16)

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

(17)

(18)

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

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

1.3.1 水平集方程的近似

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

(19)

式中：

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

(20)

(21)

(22)

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

(23)

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

(24)

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

(25)

1.3.2 重新初始化方程的近似

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

(26)

(27)

2 反演算例

2.1 反演流程

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

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

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

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

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

2.2 兩個(gè)水平集模型

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

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

(28)

φ2=1-

(29)

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

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

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

將本文反演結(jié)果與Li等[24]的反演結(jié)果做對(duì)比，見圖3?？梢钥闯?，本文采用基于任意四面體單元的磁法解析解算法[26]進(jìn)行高精度正演，并引入HJ-WENO格式對(duì)水平集方程進(jìn)行更新及重新初始化，反演得到的磁異常體邊界與模型吻合更好。

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

2.3 三個(gè)水平集模型

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

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

φ1=1-

(30)

φ2=1-

(31)

φ3=1-

(32)

圖4 三個(gè)不同磁化率模型及三個(gè)水平集計(jì)算結(jié)果(a)模型示意圖； (b)加噪總磁異常圖； (c)初始猜測(cè)模型； (d)三個(gè)水平集的反演模型；(e)圖d正演磁異常圖； (f)圖e與圖b的差值圖a中磁性體①、②、③的長(zhǎng)方體中心坐標(biāo)分別為(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細(xì)節(jié)展示(a)y=200m截面； (b)y=500m截面； (c)y=800m截面； (d)z=100m平面圖虛線為真實(shí)模型的輪廓線，實(shí)線為本文反演結(jié)果

2.4 深部異常體模型

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

模型包括兩個(gè)磁性長(zhǎng)方體，其中心坐標(biāo)分別為 (2000m,1500m,900m)、(2950m,2550m,825m)，大小分別為1000m×500m×400m、700m×300m×350m,極化率已知，分別為0.16、0.30SI?；诩釉肟偞女惓?shù)據(jù)(圖6b)，反演中引入兩個(gè)水平集函數(shù)φ1和φ2，其初始猜測(cè)(圖6c)分別為

φ1=1-

(33)

φ2=1-

(34)

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

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

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

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

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

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

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

φ=1-

(35)

圖8b為引入單水平集的反演結(jié)果?？梢娂词钩跏寄Ｐ椭患僭O(shè)了一個(gè)磁源體，反演結(jié)果也能準(zhǔn)確地反映有三個(gè)獨(dú)立的磁源體，說(shuō)明水平集的個(gè)數(shù)對(duì)反演結(jié)果影響不大。

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

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

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

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

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

3 結(jié)論

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

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

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

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