多參數(shù)約束磁性體三維形態(tài)反演

李金朋 范紅波 劉 利 張英堂 李志寧 武炳陽

(①93114部隊(duì)，北京100195；②陸軍工程大學(xué)石家莊校區(qū)，河北石家莊050003；③北京市遙感信息研究所，北京 100192)

0 引言

磁測(cè)數(shù)據(jù)反演主要應(yīng)用于軍事偵察(未爆彈、潛艇和水雷等)、水文、能源及工程地質(zhì)勘探等領(lǐng)域[1-4]。磁性目標(biāo)體的三維反演是利用觀測(cè)面上磁測(cè)數(shù)據(jù)，對(duì)其空間形狀、位置、磁化強(qiáng)度及磁化方向等參數(shù)進(jìn)行求解的過程，為下一步地質(zhì)解釋提供可靠的數(shù)據(jù)支持。

磁性體反演主要包括物性反演和形態(tài)反演兩方面[5]。物性反演的主要思路是將觀測(cè)面對(duì)應(yīng)的地下區(qū)域離散化為規(guī)則的長(zhǎng)方體網(wǎng)格，通過反演獲得這些網(wǎng)格的磁性參數(shù)，并最終獲得目標(biāo)體的磁性參數(shù)空間分布[6-8]。物性反演過程中，由于目標(biāo)函數(shù)固有的欠定性特征，計(jì)算結(jié)果存在多解性；同時(shí)，由于需要進(jìn)行多次模型正演計(jì)算，核函數(shù)矩陣會(huì)大大增加計(jì)算時(shí)間，影響計(jì)算效率。形態(tài)反演的主要思路是利用一定的計(jì)算準(zhǔn)則、基于觀測(cè)數(shù)據(jù)對(duì)地下多面體進(jìn)行擬合，利用多面體的形態(tài)對(duì)待測(cè)模型的三維分布進(jìn)行求解[9]。Uieda等[10]利用L2范數(shù)定義數(shù)據(jù)泛函約束函數(shù)對(duì)待測(cè)目標(biāo)進(jìn)行形態(tài)反演；曹書錦等[11]基于L1范數(shù)利用種子反演方法，提高了重力梯度張量反演精度；Uieda等[12]利用重力梯度張量，對(duì)不同待測(cè)目標(biāo)的“種子”分配不同密度值，實(shí)現(xiàn)特定目標(biāo)的反演，有效排除了非目標(biāo)場(chǎng)源的干擾。常規(guī)的形態(tài)反演方法能夠克服物性反演方法在迭代過程中對(duì)核函數(shù)矩陣進(jìn)行多次計(jì)算以及計(jì)算結(jié)果存在多解性的問題，但前提是需要獲得待測(cè)目標(biāo)的先驗(yàn)信息。然而，實(shí)際應(yīng)用中常利用經(jīng)驗(yàn)數(shù)據(jù)對(duì)模型的磁性參數(shù)和幾何參數(shù)進(jìn)行賦值，導(dǎo)致計(jì)算結(jié)果的精度下降。

針對(duì)形態(tài)反演方法中存在的問題，本文提出了基于多參數(shù)約束的磁性目標(biāo)三維形態(tài)反演方法：首先計(jì)算磁性目標(biāo)的平面位置、深度、磁化方向及磁化強(qiáng)度；然后，聯(lián)合垂直磁場(chǎng)數(shù)據(jù)Bz和磁張量數(shù)據(jù)，對(duì)磁性目標(biāo)進(jìn)行三維反演，在反演過程中，對(duì)磁性目標(biāo)待反演空間進(jìn)行劃分，建立待增長(zhǎng)區(qū)域，并計(jì)算最優(yōu)增長(zhǎng)模塊，實(shí)現(xiàn)磁性目標(biāo)的三維重建。

1 方法原理

本文計(jì)算過程如圖1所示。首先，根據(jù)觀測(cè)面的磁測(cè)數(shù)據(jù)，分別對(duì)磁性目標(biāo)的水平位置、深度、磁化方向和磁化強(qiáng)度進(jìn)行估計(jì)，確定初始種子的磁性參數(shù)和幾何參數(shù)；然后，利用形態(tài)反演迭代計(jì)算方法，通過確定待增長(zhǎng)區(qū)域中的最優(yōu)增長(zhǎng)模塊，實(shí)現(xiàn)模型反演；最終，將滿足終止條件的計(jì)算結(jié)果輸出，即最終反演結(jié)果。

圖1 本文算法流程圖

1.1 形態(tài)反演理論

假設(shè)觀測(cè)面為平面，d為觀測(cè)面上的實(shí)際觀測(cè)數(shù)據(jù)，b為預(yù)測(cè)模型在觀測(cè)面上的正演數(shù)據(jù)。若數(shù)據(jù)包含較大誤差，L2范數(shù)比L1范數(shù)更敏感[11，13]，因此，利用L1范數(shù)能夠獲得更加穩(wěn)定的反演結(jié)果。基于L1范數(shù)定義的數(shù)據(jù)泛函約束函數(shù)為

(1)

式中：i代表觀測(cè)點(diǎn)序號(hào)；N代表觀測(cè)面上總觀測(cè)點(diǎn)數(shù)；di、bi分別為d、b的元素。磁梯度張量矩陣可表示為

(2)

式中：Bx、By和Bz分別為磁場(chǎng)的x、y和z分量；Bαβ(α，β=x，y，z)代表磁梯度張量B的分量。

假設(shè)反演過程中需同時(shí)對(duì)H類數(shù)據(jù)進(jìn)行計(jì)算，則總數(shù)據(jù)約束函數(shù)Ψ(M)可以表示為

(3)

式中：M為反演模型的磁化強(qiáng)度；φl為第l類數(shù)據(jù)類型的約束函數(shù)。例如，當(dāng)利用磁梯度張分量Bxz、Byz、Bzz進(jìn)行反演時(shí)，H=3，φ1(M)=φxz(M)，φ2(M)=φyz(M)，φ3(M)=φzz(M)，可利用式(1)分別對(duì)其進(jìn)行計(jì)算。

形態(tài)反演過程中，包含如下約束條件[14-15]：①模型的解是連續(xù)的(內(nèi)部沒有空洞)；②各網(wǎng)格的磁化強(qiáng)度只能為M或者0；③模型增長(zhǎng)過程中，需要設(shè)置初始位置，且新增長(zhǎng)模塊與當(dāng)前模塊至少有一個(gè)平面是共面的。根據(jù)上述約束條件，形態(tài)反演目標(biāo)函數(shù)Γ(M)可以表示為

Γ(M)=Ψ(M)+ρθ(M)

(4)

式中：θ(M)為模型在空間上定義的正則化約束函數(shù)；ρ為正則化參數(shù)，用于調(diào)節(jié)Ψ(M)與θ(M)兩個(gè)函數(shù)之間的平衡。對(duì)正則化參數(shù)ρ的選擇，本文采用自適應(yīng)正則化方法

ρ(n+1)=qρ(n)

(5)

式中：n表示迭代次數(shù)；q=0.95。本文設(shè)定初始值ρ(0)=0.1。θ(M)可以表示為

(6)

式中：w、f、g分別為當(dāng)前模型在x、y、z三個(gè)方向上的長(zhǎng)度；P為當(dāng)前模型包含的長(zhǎng)方體的個(gè)數(shù)；Mj為第j個(gè)長(zhǎng)方體的磁化強(qiáng)度；ε為一個(gè)很小的正數(shù)，用于避免Mj為0時(shí)計(jì)算的不穩(wěn)定，本文定義ε=1×10-5；Lj為待增長(zhǎng)長(zhǎng)方體的中心點(diǎn)與當(dāng)前模型的第j個(gè)長(zhǎng)方體的中心點(diǎn)之間的距離。θ(M)在計(jì)算過程中加入了模型的位置屬性，通過長(zhǎng)方體單元間的距離對(duì)反演模型施加約束，避免反演結(jié)果沿某一方向無限增長(zhǎng)，可提高反演結(jié)果的準(zhǔn)確性[16]。

形態(tài)反演的具體計(jì)算過程如圖2所示。首先，在待測(cè)區(qū)域內(nèi)設(shè)置初始模型作為種子，即分別對(duì)初始模型的位置、磁化方向以及磁化強(qiáng)度進(jìn)行賦值；然后，利用式(4)計(jì)算待增長(zhǎng)區(qū)域內(nèi)所有網(wǎng)格的目標(biāo)函數(shù)值，將待增長(zhǎng)區(qū)域內(nèi)貢獻(xiàn)最大的網(wǎng)格作為增長(zhǎng)網(wǎng)格(即確保Ψ(M)變小的前提下Γ(M)最小)，并設(shè)置增長(zhǎng)模型的磁化強(qiáng)度與初始模型一致；最后，反復(fù)迭代，直到目標(biāo)函數(shù)Γ(M)達(dá)到閾值或者達(dá)到最大迭代次數(shù)，終止迭代，獲得最終反演結(jié)果。

圖2 形態(tài)反演模型增長(zhǎng)示意圖

1.2 磁性目標(biāo)多參數(shù)反演

在形態(tài)反演過程中，獲得磁性目標(biāo)的準(zhǔn)確初始磁性參數(shù)和幾何參數(shù)，對(duì)計(jì)算結(jié)果有較大的影響。傳統(tǒng)的形態(tài)反演方法基于先驗(yàn)信息對(duì)初始磁性參數(shù)和幾何參數(shù)進(jìn)行賦值，這會(huì)導(dǎo)致對(duì)先驗(yàn)信息的依賴性較強(qiáng)，對(duì)于部分先驗(yàn)信息未知的目標(biāo)，計(jì)算結(jié)果精度低。為了提高方法的適應(yīng)性，本文采用下述方案分別對(duì)磁性目標(biāo)的水平位置、深度、磁化方向以及磁化強(qiáng)度進(jìn)行估計(jì)。

1.2.1 平面位置估計(jì)

歸一化磁源強(qiáng)度(NSS)是基于磁偶極子磁梯度張量矩陣計(jì)算得到的，在剩余磁化條件下，能夠?qū)Υ判阅繕?biāo)的水平位置進(jìn)行有效計(jì)算。NSS的最大值點(diǎn)即對(duì)應(yīng)磁偶極子的實(shí)際位置。對(duì)于多邊形磁性體而言(如水平薄板)，當(dāng)測(cè)量面距離多邊形較近時(shí)(即構(gòu)造指數(shù)為小于3的正整數(shù)，模型不能等效為磁偶極子)，NSS極大值點(diǎn)即對(duì)應(yīng)磁性目標(biāo)的邊界(存在多個(gè)極大值點(diǎn))[17]。因此，需對(duì)觀測(cè)面磁測(cè)數(shù)據(jù)向上延拓，當(dāng)延拓面距離模型實(shí)際位置足夠大時(shí)，模型可以等效為磁偶極子。此時(shí)，NSS的極大值唯一，且與磁性目標(biāo)的中心位置一致。對(duì)于多個(gè)磁性目標(biāo)而言，需要根據(jù)NSS的分布特征劃分計(jì)算區(qū)域，使每個(gè)計(jì)算區(qū)域內(nèi)僅包含一個(gè)磁性目標(biāo)，然后再利用上述方法分別進(jìn)行延拓，并最終獲得所有磁性目標(biāo)的中心位置。

(7)

因此，本文基于NSS數(shù)據(jù)的特性，將觀測(cè)面上各計(jì)算區(qū)域內(nèi)的NSS極大值點(diǎn)作為磁性目標(biāo)的水平中心位置。當(dāng)計(jì)算區(qū)域內(nèi)包含多個(gè)極值時(shí)，利用向上延拓理論，將觀測(cè)數(shù)據(jù)向上延拓[20]，直到計(jì)算區(qū)域內(nèi)僅包含一個(gè)極大值點(diǎn)。

1.2.2 深度估計(jì)

謝汝寬等[21]提出了最小反演擬合差的深度估計(jì)方法，利用空間域單層等效源估計(jì)磁性目標(biāo)的深度。本文在此基礎(chǔ)上，提出了頻率域單層等效源深度估計(jì)方法，可進(jìn)一步節(jié)省計(jì)算時(shí)間。

假設(shè)觀測(cè)面為水平面，將觀測(cè)面以下空間劃分為不同深度、相同厚度的水平長(zhǎng)方體層，且長(zhǎng)方體層尺寸相同。假設(shè)某一長(zhǎng)方體層的頂面深度為z1，底面深度為z2，觀測(cè)面的高度設(shè)置為z0，此長(zhǎng)方體層的磁異常分量Bz與磁化強(qiáng)度M在頻率域的關(guān)系為[22]

(8)

以第c層為例，假設(shè)其頂面深度為zc，其磁化強(qiáng)度反演結(jié)果Mc可表示為[23]

(9)

式中h(kx，ky，zc)=(zce-kzc)s，正整數(shù)s∈[1，10]，取值越大成像結(jié)果的分辨率越高，本文設(shè)定s=10。

基于最小反演擬合差的頻率域單層等效源深度估計(jì)方法的計(jì)算過程[21]描述如下：在觀測(cè)平面下設(shè)置某一厚度的長(zhǎng)方體等效源層，該長(zhǎng)方體層的厚度為z2-z1；將等效源層向下移動(dòng)一定的距離，分別利用式(8)和式(9)計(jì)算不同等效源層在觀測(cè)平面的正演數(shù)據(jù)與實(shí)際數(shù)據(jù)的反演擬合差；在不斷向下移動(dòng)過程中，選擇擬合差最小時(shí)對(duì)應(yīng)的長(zhǎng)方體層深度作為模型的中心深度。對(duì)于多目標(biāo)體，將觀測(cè)面上的磁測(cè)數(shù)據(jù)劃分為多個(gè)僅包含一個(gè)目標(biāo)的計(jì)算區(qū)域，分別對(duì)不同計(jì)算區(qū)域內(nèi)的待測(cè)模型利用本文深度計(jì)算方法進(jìn)行計(jì)算，實(shí)現(xiàn)對(duì)不同深度(特殊情形下也可能是同一深度)的多個(gè)磁性目標(biāo)體的深度逐一進(jìn)行估計(jì)。

1.2.3 磁化方向估計(jì)

在磁化方向未知的情況下，對(duì)磁性目標(biāo)的磁傾角和磁偏角等間隔選取一系列的數(shù)據(jù)點(diǎn)，組成磁化方向暫定值。對(duì)這一系列磁化方向下的化極異常值與歸一化磁源強(qiáng)度進(jìn)行試錯(cuò)，得到不同的互相關(guān)系數(shù)，當(dāng)互相關(guān)系數(shù)取得最大值時(shí)，即表明對(duì)應(yīng)的磁化方向?yàn)樽罴压烙?jì)值。

定義實(shí)測(cè)區(qū)域內(nèi)磁異常的化極結(jié)果ΔTrtp與歸一化磁源強(qiáng)度NSS的互相關(guān)系數(shù)為

(10)

1.2.4 磁化強(qiáng)度估計(jì)

在傳統(tǒng)磁性目標(biāo)形態(tài)反演方法中，初始磁化強(qiáng)度根據(jù)先驗(yàn)信息確定，為了更準(zhǔn)確地獲得磁性目標(biāo)的磁化強(qiáng)度，本文提出下述磁化強(qiáng)度估計(jì)方法。

(1)估計(jì)磁性目標(biāo)初始種子的位置及磁化方向信息，給定初始磁化強(qiáng)度，并對(duì)磁性目標(biāo)進(jìn)行反演。

(2)根據(jù)步驟(1)的反演結(jié)果，結(jié)合磁性目標(biāo)的位置信息及磁化方向信息，將磁化強(qiáng)度按照一定的間隔，取一系列值，并按照從小到大的順序進(jìn)行正演計(jì)算。

(3)定義均方根誤差

計(jì)算不同磁化強(qiáng)度條件下正演數(shù)據(jù)與觀測(cè)數(shù)據(jù)的RMSE，并將RMSE最小值對(duì)應(yīng)的磁化強(qiáng)度作為當(dāng)前模型的磁化強(qiáng)度；

(4)重復(fù)步驟(1)～步驟(3)，并將初始磁化強(qiáng)度定義為上一次計(jì)算得到的最優(yōu)磁化強(qiáng)度值，直到RMSE小于預(yù)設(shè)的精度或者達(dá)到預(yù)設(shè)的迭代次數(shù)。

1.3 形態(tài)反演迭代終止準(zhǔn)則

在形態(tài)反演過程中，增長(zhǎng)模型為目標(biāo)函數(shù)Γ(M)最小時(shí)對(duì)應(yīng)的長(zhǎng)方體模型。定義

(11)

式中Φold(M)、Φnew(M)分別為未加入和已加入最優(yōu)增長(zhǎng)模型的數(shù)據(jù)約束函數(shù)。通過計(jì)算η(M)判斷是否終止迭代，當(dāng)η(M)小于預(yù)設(shè)的計(jì)算精度時(shí)則停止計(jì)算，這個(gè)值一般設(shè)定為1×10-4～1×10-6。

2 仿真分析

2.1 長(zhǎng)方體模型

為了證明本文方法的有效性，建立圖3所示的長(zhǎng)方體模型進(jìn)行仿真分析。假設(shè)觀測(cè)面的深度為0(即地面)，觀測(cè)點(diǎn)距為2m，觀測(cè)區(qū)域?yàn)?6m×26m，中心點(diǎn)與坐標(biāo)原點(diǎn)重合。長(zhǎng)方體模型的長(zhǎng)、寬、高分別為10、10、6m，中心位置坐標(biāo)為(0，0，12m)。磁性體的磁化強(qiáng)度為40A/m，磁傾角I=70°，磁偏角D=20°，背景磁化方向與感應(yīng)磁化方向相同。

加入信噪比為40%的高斯噪聲，該模型觀測(cè)面上的磁測(cè)數(shù)據(jù)如圖4所示。利用本文方法對(duì)長(zhǎng)方體的磁性參數(shù)和幾何參數(shù)進(jìn)行反演，結(jié)果如圖5所示。由圖5可知，計(jì)算得到的磁性體中心位置為(1m，1m，12m)，磁傾角I=65°，磁偏角D=19°。將此計(jì)算結(jié)果作為反演的初始種子參數(shù)值。按照前文敘述，在不同觀測(cè)數(shù)據(jù)組合條件下對(duì)磁性目標(biāo)進(jìn)行反演。設(shè)磁化強(qiáng)度為40A/m，分別利用不同的數(shù)據(jù)組合，即Bzz、Bxx/Bxy/Bxz/Byy/Byz/Bzz、Bz/Bzz以及Bz/Bxx/Bxy/Bxz/Byy/Byz/Bzz進(jìn)行反演，反演結(jié)果如圖6所示。

對(duì)比利用四種不同的組合數(shù)據(jù)反演結(jié)果(圖6)，可以看出，相較于加入Bz數(shù)據(jù)的反演結(jié)果(圖6c、圖6d)，單獨(dú)利用Bzz數(shù)據(jù)(圖a)或Bxx/Bxy/Bxz/Byy/Byz/Bzz(圖b)數(shù)據(jù)組合獲得的反演結(jié)果，其垂直分辨率較低。

圖3 長(zhǎng)方體模型示意圖

圖4 長(zhǎng)方體模型加噪磁張量及NSS計(jì)算結(jié)果(a)Bxx； (b) Bxy； (c)Bxz； (d)Byy； (e)Byz； (f)Bzz； (g)Bz； (h)NSS(白線方框?yàn)榇判泽w的水平位置)

圖5 長(zhǎng)方體模型參數(shù)估計(jì)結(jié)果(a)磁化方向； (b)水平位置； (c)深度

對(duì)不同數(shù)據(jù)獲得的反演模型進(jìn)行正演計(jì)算，計(jì)算結(jié)果的RMSE見表1。可以看出，若反演數(shù)據(jù)組合中包含Bz，反演模型的正演張量數(shù)據(jù)精度稍有下降。根據(jù)圖6和表1可以看出，雖然加入Bz數(shù)據(jù)后反演模型的正演數(shù)據(jù)精度稍有下降，但是利用Bz/Bzz數(shù)據(jù)組合和Bz/Bxx/Bxy/Bxz/Byy/Byz/Bzz數(shù)據(jù)組合能夠獲得更加準(zhǔn)確的三維反演結(jié)果。同時(shí)，與Bxx/Bxy/Bxz/Byy/Byz/Bzz數(shù)據(jù)組合相比，利用Bz/Bzz數(shù)據(jù)組合的反演效率更高。

利用獲得的中心位置作為種子的初始位置，分別假設(shè)初始磁化強(qiáng)度為1、200A/m，利用Bz/Bzz進(jìn)行反演，得到不同初始磁化強(qiáng)度下模型的反演結(jié)果(圖7a)。利用本文的磁化強(qiáng)度估計(jì)方法計(jì)算圖7a模型的磁化強(qiáng)度，獲得的磁化強(qiáng)度結(jié)果如圖7b所示。可以看出，即便設(shè)置不同的磁化強(qiáng)度作為初始種子，反演獲得的磁化強(qiáng)度與實(shí)際值均比較接近，說明初始磁化強(qiáng)度的大小對(duì)反演結(jié)果影響不大。

圖6 不同數(shù)據(jù)組合條件下的磁性體形態(tài)反演結(jié)果

表1 不同數(shù)據(jù)組合條件下反演結(jié)果的正演數(shù)據(jù)RMSE統(tǒng)計(jì)

圖7 初始磁化強(qiáng)度為1A/m(左)和200A/m(右)條件下模型三維反演結(jié)果(a)及設(shè)定不同初始磁化強(qiáng)度下反演模型的正演誤差(b)

2.2 傾斜板狀體

為了分析估計(jì)的位置參數(shù)對(duì)反演結(jié)果的影響，分別設(shè)定三個(gè)不同初始位置，利用Bz/Bzz數(shù)據(jù)進(jìn)行反演計(jì)算。初始位置1位于上頂面邊緣的中心位置，初始位置2即本文方法獲得的位置，初始位置3位于下底面邊緣的中心位置(圖11上)。圖11(下)為對(duì)應(yīng)的反演結(jié)果。可以看出，本文方法在不同的初始位置條件下都能夠?qū)Υ判阅繕?biāo)的形態(tài)進(jìn)行有效反演。但是，當(dāng)初始位置設(shè)定在上頂面邊緣(圖11a)時(shí)，反演結(jié)果更集中于上頂層部分；而初始位置位于下底面邊緣(圖11c)時(shí)，反演結(jié)果更集中于底層。分別對(duì)這三種反演結(jié)果進(jìn)行正演計(jì)算，計(jì)算結(jié)果的均方根誤差RMSE如表2所示。根據(jù)表2可以看出，相較于NSS和TMI，以本文方法獲得的位置(圖11b)作為模型初始種子位置進(jìn)行反演時(shí)，結(jié)果具有更高的計(jì)算精度。

圖8 傾斜板狀磁性體模型空間示意圖

進(jìn)一步地，對(duì)NSS數(shù)據(jù)和TMI數(shù)據(jù)的反演能力進(jìn)行討論。分別利用NSS數(shù)據(jù)和TMI數(shù)據(jù)進(jìn)行形態(tài)反演，結(jié)果如圖12所示。對(duì)反演模型進(jìn)行正演計(jì)算，其結(jié)果與原始模型正演數(shù)據(jù)RMSE統(tǒng)計(jì)如表2所示。由圖12可知，相較于TMI，NSS反演結(jié)果在深度方向上的分辨率較低，這是由于NSS數(shù)據(jù)與距離呈4次方衰減，深度信息衰減較快，而TMI數(shù)據(jù)與距離呈3次方衰減。此外，相較于Bz/Bzz組合數(shù)據(jù)，由于NSS和TMI數(shù)據(jù)僅包含振幅信息，反演模型的分辨率較低。

圖9 傾斜板狀磁性體模型正演磁測(cè)數(shù)據(jù)(加噪)(a)Bz； (b)Bzz； (c)TMI； (d)NSS。圖中白線為磁性目標(biāo)水平位置

圖10 傾斜板狀磁性體模型本文方法參數(shù)估計(jì)結(jié)果(a) 磁化方向； (b)水平位置； (c)深度； (d) 磁化強(qiáng)度

圖11 傾斜板狀磁性體模型種子點(diǎn)初始位置1(a)、2(b)、3(c)的反演結(jié)果上圖為種子點(diǎn)位置，下圖為對(duì)應(yīng)的反演結(jié)果

表2 不同初始位置條件下不同數(shù)據(jù)組合反演模型的不同參數(shù)正演結(jié)果RMSE統(tǒng)計(jì)

圖12 傾斜板狀磁性體模型NSS(上)和TMI(下)反演結(jié)果

2.3 多磁性體模型

在實(shí)際計(jì)算過程中，模型會(huì)受剩余磁化的影響，導(dǎo)致實(shí)際磁化方向與背景磁化方向不同。因此，若直接利用背景場(chǎng)的磁化方向進(jìn)行計(jì)算，會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生一定的影響，進(jìn)而影響反演結(jié)果[3]。為了解決這一問題，常使用弱敏感于磁化方向的NSS以及TMI進(jìn)行反演。

建立一個(gè)傾斜板狀體和長(zhǎng)方體的組合模型(圖13，表3)，對(duì)剩磁條件下不同磁測(cè)數(shù)據(jù)的反演能力進(jìn)行分析。背景磁化方向?yàn)镮=60°，D=20°。觀測(cè)面深度為0(即位于地面)，觀測(cè)點(diǎn)間距為2m，觀測(cè)區(qū)域大小為26m×26m。觀測(cè)數(shù)據(jù)中加入40%的高斯噪聲，觀測(cè)面上的正演磁測(cè)數(shù)據(jù)如圖14所示。

基于NSS主正極值分布范圍與磁性目標(biāo)的垂直分布范圍基本一致，且受剩余磁化影響較小的特點(diǎn)，分別研究分析長(zhǎng)方體目標(biāo)和傾斜板狀目標(biāo)，其對(duì)應(yīng)的區(qū)域如圖14d中黑線所示。利用本文的磁性目標(biāo)多參數(shù)反演方法對(duì)模型初始參數(shù)進(jìn)行計(jì)算，獲得的初始位置、磁化方向和磁化強(qiáng)度如圖15所示。根據(jù)圖15上可知，計(jì)算得到的長(zhǎng)方體模型中心位置坐標(biāo)為(1m，-10m，10m)，磁化方向?yàn)镮=74°，D=21°，磁化強(qiáng)度為38A/m；根據(jù)圖15下，傾斜板模型的中心位置坐標(biāo)為(1m，13m，8m)，磁化方向?yàn)镮=44°，D=31°，磁化強(qiáng)度為22A/m。對(duì)比實(shí)際模型參數(shù)可知，對(duì)于多目標(biāo)體模型的參數(shù)估計(jì)，由于各磁性體的相互影響，模型的參數(shù)估計(jì)精度有所下降。

圖13 組合模型示意圖

圖14 組合模型磁測(cè)數(shù)據(jù)

表3 多磁性體模型參數(shù)

圖15 利用本文方法得到的區(qū)域1(上)和區(qū)域2(下)的磁化方向(a)、水平位置(b)、深度(c)和磁化強(qiáng)度(d)估計(jì)結(jié)果

圖16 剩磁條件下利用Bz/Bzz(a)、TMI(b)和NSS(c)反演的多目標(biāo)體模型空間形態(tài)

利用本文方法估計(jì)組合模型的初始位置、磁化方向及磁化強(qiáng)度(圖15)，分別利用Bz/Bzz、TMI以及NSS對(duì)此模型的空間位置進(jìn)行反演，結(jié)果如圖16所示。可以看出，在剩余磁化條件下，與TMI和NSS相比，基于Bz/Bzz數(shù)據(jù)的反演結(jié)果受相鄰磁性目標(biāo)的影響最小，具有更高的水平及垂向反演精度。

3 試驗(yàn)驗(yàn)證

在河北省石家莊市某地對(duì)一圓柱狀磁性體目標(biāo)進(jìn)行實(shí)測(cè)，以驗(yàn)證本文方法。裝置系統(tǒng)(圖17)主要包括Mag-03傳感器、數(shù)字采集模塊及軟件操作終端。實(shí)驗(yàn)中，傳感器固定在無磁實(shí)驗(yàn)臺(tái)架上，采用掃描方式對(duì)待測(cè)區(qū)域內(nèi)的每一測(cè)點(diǎn)進(jìn)行測(cè)量。測(cè)區(qū)大小為1.9m×1.9m，測(cè)點(diǎn)間隔為0.1m。在測(cè)區(qū)內(nèi)放置一個(gè)南北走向的水平圓柱體鐵塊(磁性體)。背景磁化方向?yàn)镮=56°，D=-16°；圓柱體底面直徑為0.10m，軸長(zhǎng)為1.05m，中心位置坐標(biāo)為(0.80m，0.95m，0.45m)。觀測(cè)面的實(shí)際磁測(cè)數(shù)據(jù)如圖18所示。

圖17 試驗(yàn)裝置及待測(cè)模型

圖18 實(shí)驗(yàn)實(shí)測(cè)數(shù)據(jù)(a)Bxx； (b) Bxy； (c)Bxz； (d)Byy； (e) Byz； (f)Bzz； (g)Bz； (h)NSS； (i)TMI

利用本文方法對(duì)模型初始參數(shù)進(jìn)行計(jì)算，獲得的異常體初始位置、磁化方向和磁化強(qiáng)度結(jié)果如圖19所示。基于實(shí)測(cè)數(shù)據(jù)計(jì)算得到的模型中心點(diǎn)坐標(biāo)為(0.8m，0.9m，0.4m)，磁化方向?yàn)镮=21°，D=-1°，磁化強(qiáng)度為960A/m。利用這些模型磁性參數(shù)和幾何參數(shù)，分別基于Bz/Bzz、TMI以及NSS對(duì)待測(cè)模型進(jìn)行反演，結(jié)果如圖20所示。再對(duì)這些反演模型進(jìn)行正演，其RMSE數(shù)據(jù)統(tǒng)計(jì)如表4所示。根據(jù)表4及圖20可以看出，基于Bz/Bzz數(shù)據(jù)獲得的計(jì)算結(jié)果與模型實(shí)際位置最接近，具有較高的計(jì)算精度。

表4 不同數(shù)據(jù)反演模型的正演數(shù)據(jù)RMSE統(tǒng)計(jì)

圖19 基于實(shí)測(cè)數(shù)據(jù)的模型參數(shù)估計(jì)結(jié)果(a)磁化方向； (b)水平位置； (c)深度； (d)磁化強(qiáng)度

圖20 基于Bz/Bzz(a)、TMI(b)和NSS(c)反演的磁性體空間形態(tài)上圖為鳥瞰圖，中圖為側(cè)視圖，下圖為三維顯示

4 結(jié)論

本文針對(duì)磁性目標(biāo)三維形態(tài)反演中存在的問題，提出了一種磁性目標(biāo)三維形態(tài)反演方法，即根據(jù)實(shí)測(cè)磁數(shù)據(jù)，對(duì)磁性目標(biāo)的水平位置、深度、磁化方向和磁化強(qiáng)度進(jìn)行初步估計(jì)，確定初始種子的磁性參數(shù)和幾何參數(shù)；然后利用形態(tài)反演迭代計(jì)算方法，通過確定待增長(zhǎng)區(qū)域中的最優(yōu)增長(zhǎng)模塊，反演得到模型空間分布形態(tài)。仿真和實(shí)驗(yàn)結(jié)果證明了本文方法的正確性和實(shí)用性。具體得到如下結(jié)論。

(1)對(duì)于磁性體的形態(tài)反演，聯(lián)合Bz與磁張量，能夠改善磁張量在深度上反演分辨率低的問題。

(2)本文方法無需先驗(yàn)信息，只需通過觀測(cè)數(shù)據(jù)就可以直接計(jì)算磁性目標(biāo)體的中心位置、磁化方向和磁化強(qiáng)度，提高了形態(tài)反演方法的適用性，減小了人為經(jīng)驗(yàn)確定反演初始值所引入的誤差。

(3)在剩余磁化條件下，本文方法能夠直接利用Bz/Bzz數(shù)據(jù)組合對(duì)多目標(biāo)模型進(jìn)行反演，與基于NSS和TMI的反演結(jié)果相比，明顯提高了反演分辨率。

利用NSS數(shù)據(jù)對(duì)磁性目標(biāo)的水平位置進(jìn)行估計(jì)時(shí)，臨近疊加異常會(huì)產(chǎn)生混疊，影響計(jì)算精度；對(duì)于臺(tái)階等復(fù)雜模型，向上延拓的中心可能與磁性目標(biāo)的中心位置不一致。因此，在下一步工作中需要對(duì)上述問題進(jìn)一步研究。