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

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

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

0 引言

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

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

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

1 方法原理

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

圖1 本文算法流程圖

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

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

(1)

式中：i代表觀測點序號；N代表觀測面上總觀測點數；di、bi分別為d、b的元素。磁梯度張量矩陣可表示為

(2)

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

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

(3)

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

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

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

(4)

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

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

(5)

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

(6)

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

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

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

1.2 磁性目標多參數反演

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

1.2.1 平面位置估計

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

(7)

因此，本文基于NSS數據的特性，將觀測面上各計算區(qū)域內的NSS極大值點作為磁性目標的水平中心位置。當計算區(qū)域內包含多個極值時，利用向上延拓理論，將觀測數據向上延拓[20]，直到計算區(qū)域內僅包含一個極大值點。

1.2.2 深度估計

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

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

(8)

以第c層為例，假設其頂面深度為zc，其磁化強度反演結果Mc可表示為[23]

(9)

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

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

1.2.3 磁化方向估計

在磁化方向未知的情況下，對磁性目標的磁傾角和磁偏角等間隔選取一系列的數據點，組成磁化方向暫定值。對這一系列磁化方向下的化極異常值與歸一化磁源強度進行試錯，得到不同的互相關系數，當互相關系數取得最大值時，即表明對應的磁化方向為最佳估計值。

定義實測區(qū)域內磁異常的化極結果ΔTrtp與歸一化磁源強度NSS的互相關系數為

(10)

1.2.4 磁化強度估計

在傳統磁性目標形態(tài)反演方法中，初始磁化強度根據先驗信息確定，為了更準確地獲得磁性目標的磁化強度，本文提出下述磁化強度估計方法。

(1)估計磁性目標初始種子的位置及磁化方向信息，給定初始磁化強度，并對磁性目標進行反演。

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

(3)定義均方根誤差

計算不同磁化強度條件下正演數據與觀測數據的RMSE，并將RMSE最小值對應的磁化強度作為當前模型的磁化強度；

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

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

在形態(tài)反演過程中，增長模型為目標函數Γ(M)最小時對應的長方體模型。定義

(11)

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

2 仿真分析

2.1 長方體模型

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

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

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

圖3 長方體模型示意圖

圖4 長方體模型加噪磁張量及NSS計算結果(a)Bxx； (b) Bxy； (c)Bxz； (d)Byy； (e)Byz； (f)Bzz； (g)Bz； (h)NSS(白線方框為磁性體的水平位置)

圖5 長方體模型參數估計結果(a)磁化方向； (b)水平位置； (c)深度

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

利用獲得的中心位置作為種子的初始位置，分別假設初始磁化強度為1、200A/m，利用Bz/Bzz進行反演，得到不同初始磁化強度下模型的反演結果(圖7a)。利用本文的磁化強度估計方法計算圖7a模型的磁化強度，獲得的磁化強度結果如圖7b所示?？梢钥闯?，即便設置不同的磁化強度作為初始種子，反演獲得的磁化強度與實際值均比較接近，說明初始磁化強度的大小對反演結果影響不大。

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

表1 不同數據組合條件下反演結果的正演數據RMSE統計

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

2.2 傾斜板狀體

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

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

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

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

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

圖11 傾斜板狀磁性體模型種子點初始位置1(a)、2(b)、3(c)的反演結果上圖為種子點位置，下圖為對應的反演結果

表2 不同初始位置條件下不同數據組合反演模型的不同參數正演結果RMSE統計

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

2.3 多磁性體模型

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

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

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

圖13 組合模型示意圖

圖14 組合模型磁測數據

表3 多磁性體模型參數

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

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

利用本文方法估計組合模型的初始位置、磁化方向及磁化強度(圖15)，分別利用Bz/Bzz、TMI以及NSS對此模型的空間位置進行反演，結果如圖16所示?？梢钥闯觯谑Ｓ啻呕瘲l件下，與TMI和NSS相比，基于Bz/Bzz數據的反演結果受相鄰磁性目標的影響最小，具有更高的水平及垂向反演精度。

3 試驗驗證

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

圖17 試驗裝置及待測模型

圖18 實驗實測數據(a)Bxx； (b) Bxy； (c)Bxz； (d)Byy； (e) Byz； (f)Bzz； (g)Bz； (h)NSS； (i)TMI

利用本文方法對模型初始參數進行計算，獲得的異常體初始位置、磁化方向和磁化強度結果如圖19所示?；趯崪y數據計算得到的模型中心點坐標為(0.8m，0.9m，0.4m)，磁化方向為I=21°，D=-1°，磁化強度為960A/m。利用這些模型磁性參數和幾何參數，分別基于Bz/Bzz、TMI以及NSS對待測模型進行反演，結果如圖20所示。再對這些反演模型進行正演，其RMSE數據統計如表4所示。根據表4及圖20可以看出，基于Bz/Bzz數據獲得的計算結果與模型實際位置最接近，具有較高的計算精度。

表4 不同數據反演模型的正演數據RMSE統計

圖19 基于實測數據的模型參數估計結果(a)磁化方向； (b)水平位置； (c)深度； (d)磁化強度

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

4 結論

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

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

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

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

利用NSS數據對磁性目標的水平位置進行估計時，臨近疊加異常會產生混疊，影響計算精度；對于臺階等復雜模型，向上延拓的中心可能與磁性目標的中心位置不一致。因此，在下一步工作中需要對上述問題進一步研究。