李金朋， 張英堂， 范紅波， 李志寧， 張光

(1.陸軍工程大學(xué)石家莊校區(qū) 車(chē)輛與電氣工程系， 河北 石家莊 050003； 2.65185部隊(duì)， 遼寧 鐵嶺 112611)

0 引言

磁法勘探主要應(yīng)用于軍事偵察(未爆彈、潛艇和水雷等)、水文及工程地質(zhì)勘探等領(lǐng)域[1-2]。反演是解釋磁性目標(biāo)的重要一步，其主要目的是利用獲得的磁性目標(biāo)磁測(cè)數(shù)據(jù)來(lái)評(píng)估地下或水下未知磁性體的輪廓形態(tài)[3-4]。在反演過(guò)程中，常常假設(shè)磁性體僅包含感應(yīng)磁化而忽略剩余磁化影響。然而，在實(shí)際測(cè)量中，磁性體除了會(huì)受到感應(yīng)磁化，還會(huì)受到剩余磁化的影響，導(dǎo)致總磁化方向產(chǎn)生改變，進(jìn)而影響磁測(cè)數(shù)據(jù)的解釋精度[5-6]。

針對(duì)強(qiáng)剩磁條件下的反演方法，國(guó)內(nèi)外學(xué)者進(jìn)行了一系列相關(guān)研究。第一種方法是首先估計(jì)磁化方向，然后利用磁化方向估計(jì)值對(duì)磁性體進(jìn)行反演[7-8]。該方法的缺點(diǎn)是僅適用于孤立磁源，對(duì)于多目標(biāo)磁源則無(wú)法利用此方法進(jìn)行反演。第二種方法是將獲得的磁異常值轉(zhuǎn)化為受剩余磁化影響小、與場(chǎng)源位置對(duì)應(yīng)良好的轉(zhuǎn)換量，然后利用轉(zhuǎn)換量進(jìn)行后續(xù)反演[9-11]。上述基于轉(zhuǎn)換量的反演方法不僅適用于孤立場(chǎng)源，也適用于磁化方向不同的多個(gè)場(chǎng)源體，缺點(diǎn)是由于缺乏相位信息，僅能反演出場(chǎng)源體的大致位置。同時(shí)，在反演過(guò)程中，由于實(shí)際測(cè)量數(shù)據(jù)中存在噪聲和固有的場(chǎng)源不唯一性，導(dǎo)致磁性目標(biāo)的反演是病態(tài)多解的。為了減少測(cè)量誤差等因素的影響，常常使用L2或Tikhonov正則化思想[12-13]。吳小平等[14]提出將共軛梯度法應(yīng)用于電阻率三維反演；朱自強(qiáng)等[15]利用混合正則化反演方法減弱了反演的聚焦效應(yīng)，能有效獲得異常邊界；楊嬌嬌等[16]通過(guò)向目標(biāo)函數(shù)中引入深度加權(quán)函數(shù)來(lái)提高對(duì)目標(biāo)深度信息的識(shí)別能力；秦朋波等[17]利用再加權(quán)的方法實(shí)現(xiàn)深度加權(quán)，對(duì)目標(biāo)進(jìn)行反演；Wang等[18]利用改進(jìn)的預(yù)處理共軛梯度法對(duì)航空重力數(shù)據(jù)進(jìn)行反演，提高了反演速度和效果。上述方法為磁測(cè)數(shù)據(jù)反演提供了較好的思路，但是在反演過(guò)程中正則化參數(shù)的選擇常采用經(jīng)驗(yàn)定值，導(dǎo)致計(jì)算精度較低，難以達(dá)到最佳的擬合效果。

本文針對(duì)上述反演過(guò)程中存在的問(wèn)題，提出了強(qiáng)剩磁條件下磁性目標(biāo)三維正則化聚焦反演方法。通過(guò)向目標(biāo)函數(shù)中引入深度加權(quán)矩陣和最小支撐矩陣，避免了反演解的不穩(wěn)定性；對(duì)目標(biāo)函數(shù)進(jìn)行奇異值分解(SVD)，并利用無(wú)偏風(fēng)險(xiǎn)估計(jì)(UPRE)準(zhǔn)則對(duì)正則化參數(shù)進(jìn)行自適應(yīng)選擇，實(shí)現(xiàn)了迭代過(guò)程的自動(dòng)化。根據(jù)磁總場(chǎng)模量數(shù)據(jù)對(duì)磁源數(shù)量進(jìn)行判斷，在反演數(shù)據(jù)的選擇上，針對(duì)孤立磁源，利用磁化方向估計(jì)方法獲得實(shí)際磁化方向，通過(guò)磁梯度張量數(shù)據(jù)進(jìn)行反演；針對(duì)多目標(biāo)磁源，采用磁總場(chǎng)模量數(shù)據(jù)進(jìn)行反演。

1 三維正則化聚焦反演理論

1.1 Tikhonov正則化反演理論

磁性目標(biāo)正演公式為

do=Am+e，

(1)

式中：核函數(shù)矩陣A的元素Aij為第i(i=1，2，3，…，p)個(gè)觀測(cè)點(diǎn)觀測(cè)的第j(j=1，2，3，…，n)個(gè)單元的響應(yīng)；do∈Rp為觀測(cè)數(shù)據(jù)矩陣；m∈Rn為待測(cè)模型參數(shù)矩陣；e∈Rp表示實(shí)測(cè)數(shù)據(jù)中的誤差。

根據(jù)(1)式可知，磁性目標(biāo)的反演過(guò)程就是利用包含噪聲的實(shí)測(cè)數(shù)據(jù)do求解模型參數(shù)m. 利用Tikhonov正則化公式，可將目標(biāo)函數(shù)矩陣轉(zhuǎn)化為如下最小二乘公式：

(2)

(3)

(4)

(5)

因?yàn)榫仃嘍為可逆矩陣，所以

(6)

(7)

m(α)=ma+D-1J(α).

(8)

(9)

m(k+1)=m(k)+(D(k+1))-1J(α(k+1)).

(10)

(11)

式中：ui和vi分別為矩陣U和V的第i列。

1.2 UPRE正則化參數(shù)選擇方法

(12)

(13)

殘差的二范數(shù)均值

(14)

因此最優(yōu)正則化參數(shù)αo表示為

(15)

在Tikhonov正則化方法中，常常利用經(jīng)驗(yàn)定值或L-curve曲線等方法對(duì)正則化參數(shù)進(jìn)行定值選擇[12-13]，在迭代過(guò)程中正則化參數(shù)不隨迭代過(guò)程而發(fā)生改變。而本文提出的UPRE準(zhǔn)則在迭代過(guò)程中不斷地優(yōu)化和求解迭代參數(shù)，實(shí)現(xiàn)了正則化參數(shù)的自適應(yīng)選擇。

2 強(qiáng)剩磁條件下反演理論

磁梯度張量數(shù)據(jù)在反演過(guò)程中具有高的模型分辨率。但是，在強(qiáng)剩磁條件下模型的實(shí)際磁化方向難以判別，進(jìn)而影響最終的反演結(jié)果。針對(duì)孤立磁源，可以通過(guò)計(jì)算其實(shí)際的磁化方向進(jìn)行反演。當(dāng)待測(cè)模型體為多目標(biāo)磁源時(shí)，在剩磁影響下，直接利用磁梯度張量進(jìn)行計(jì)算的精度將大大降低，因此，需要采用弱敏感于磁化方向的磁總場(chǎng)模量數(shù)據(jù)來(lái)進(jìn)行反演。

2.1 孤立磁源

針對(duì)孤立磁源，本文提出了利用歸一化磁源強(qiáng)度和化極(NSSRP)互相關(guān)的方法對(duì)磁化方向進(jìn)行估計(jì)。

歸一化磁源強(qiáng)度是由磁梯度張量矩陣推導(dǎo)出來(lái)的1個(gè)旋轉(zhuǎn)不變量，其分布位于場(chǎng)源中心?；谶@種性質(zhì)，在磁化方向未知情況下將鐵磁形體的磁傾角和磁偏角等間隔地選取一系列數(shù)據(jù)點(diǎn)，組成磁化方向暫定值。利用這一系列磁化方向下的化極異常值與歸一化磁源強(qiáng)度進(jìn)行試錯(cuò)，可得到不同的互相關(guān)系數(shù)，其中互相關(guān)系數(shù)最大值對(duì)應(yīng)的磁化方向即為磁化方向估計(jì)值。

定義實(shí)測(cè)區(qū)域內(nèi)的磁異?；瘶O結(jié)果ΔTr與歸一化磁源強(qiáng)度數(shù)值Tn的互相關(guān)系數(shù)C為

(16)

2.2 多目標(biāo)磁源

針對(duì)多目標(biāo)磁源，本文利用磁總場(chǎng)模量進(jìn)行反演。磁總場(chǎng)模量數(shù)據(jù)為弱敏感于磁化方向的轉(zhuǎn)換量，其定義為

(17)

式中：Ba為總場(chǎng)磁異常矢量，|Ba|為磁異常矢量的模量；Bx、By和Bz分別為磁場(chǎng)矢量在x軸、y軸、z軸方向的分量。

由(17)式可以看出，磁異常模量與地下待測(cè)模型間呈非線性關(guān)系，因此需要計(jì)算其核函數(shù)矩陣。在模型進(jìn)行第k次迭代時(shí)，Ba的第v個(gè)分量Bav對(duì)地下第w個(gè)模型參數(shù)mw的核函數(shù)矩陣Avw為

(18)

在實(shí)際應(yīng)用中，由于磁總場(chǎng)模量數(shù)據(jù)受剩余磁化方向的影響較小[7]，其主正極值與實(shí)際磁性體分布位置基本一致，根據(jù)磁性體的磁總場(chǎng)模量異常數(shù)據(jù)可以判斷出磁性體的實(shí)際分布數(shù)量。

3 數(shù)值仿真

為了便于對(duì)不同磁測(cè)數(shù)據(jù)的反演效果進(jìn)行分析，構(gòu)建一個(gè)正演模型對(duì)上述方法進(jìn)行驗(yàn)證。將地下待測(cè)空間劃分為22×22×10=4 840個(gè)單元格，每個(gè)單元格均為邊長(zhǎng)為0.1 m的正方體，地下正演模型如圖1所示。正演模型的感應(yīng)磁化方向I=60°，磁偏角D=-20°，總磁化方向I=50°，磁偏角D=135°.

3.1 磁梯度張量反演

在強(qiáng)剩磁條件下利用磁梯度張量數(shù)據(jù)進(jìn)行反演，需要已知磁性體的實(shí)際磁化方向。如圖2所示，在反演過(guò)程中向磁梯度張量分量Azz中加入0.05×std(d,1)×rand(484,1)的誤差作為實(shí)測(cè)數(shù)據(jù)，其中d表示待測(cè)數(shù)據(jù)的理論值。分別利用本文提出的NSSRP互相關(guān)系數(shù)法及Gerovska[6]提出的利用磁總場(chǎng)異常和化極(TMFRP)的互相關(guān)估計(jì)方法對(duì)磁化方向進(jìn)行估計(jì)，得到的總磁化方向估計(jì)值如表1所示。

表1 兩種方法磁化方向估計(jì)結(jié)果

從表1中可以看出，與TMFRP互相關(guān)法相比，NSSRP互相關(guān)法具有更高的求解精度，進(jìn)一步證明了本文磁化方向估計(jì)方法的有效性。利用磁化方向估計(jì)值對(duì)Azz分量進(jìn)行反演，獲得的反演結(jié)果如圖3所示，其中白色實(shí)線為異常體的實(shí)際輪廓位置。根據(jù)反演結(jié)果可知，其傾斜輪廓是清晰可見(jiàn)的，其水平和垂直分辨率較高，能夠較好地反映磁性體的實(shí)際輪廓。仿真結(jié)果表明：利用NSSRP互相關(guān)系數(shù)法能夠提高磁化方向的估計(jì)精度，利用磁梯度張量分量Azz進(jìn)行反演，獲得的結(jié)果具有較好的橫向和縱向分辨率。

3.2 磁總場(chǎng)模量反演

針對(duì)圖1所示的正演模型，根據(jù)(17)式計(jì)算正演所需要的模量數(shù)據(jù)，并利用磁總場(chǎng)模量數(shù)據(jù)進(jìn)行反演。分別計(jì)算感應(yīng)磁化條件下和總磁化條件下的磁總場(chǎng)模量數(shù)據(jù)，獲得磁總場(chǎng)模量數(shù)據(jù)如圖4所示。從圖4中可以看出，在兩種磁化方向下獲得的模量數(shù)據(jù)具有高度一致性，證明了模量數(shù)據(jù)對(duì)磁化方向的不敏感性。利用模量數(shù)據(jù)進(jìn)行反演，獲得的反演結(jié)果如圖5所示。與磁梯度張量Azz單分量反演方法相比，其橫向和縱向的模型輪廓分布較好，但是其傾斜度并不明顯。這是因?yàn)槟Ａ繑?shù)據(jù)僅包含幅值信息，并不包含相位信息。

表2 平均相對(duì)誤差LRSME

通過(guò)對(duì)兩種磁測(cè)數(shù)據(jù)的反演結(jié)果進(jìn)行分析可知：利用磁梯度張量分量Azz進(jìn)行反演獲得的結(jié)果具有較好的橫向和縱向分辨率，但是在實(shí)際應(yīng)用過(guò)程中，該方法僅能對(duì)孤立磁性目標(biāo)進(jìn)行反演；利用磁總場(chǎng)模量數(shù)據(jù)進(jìn)行反演，能夠有效規(guī)避需要已知磁化方向這一條件，在實(shí)際反演過(guò)程中更加靈活，能夠?qū)崿F(xiàn)對(duì)組合磁性目標(biāo)進(jìn)行反演，但是反演效果較Azz分量稍差。因此，在反演計(jì)算過(guò)程中，需要根據(jù)實(shí)際情況選擇合適的磁測(cè)數(shù)據(jù)進(jìn)行反演計(jì)算。

4 實(shí)驗(yàn)驗(yàn)證

為了驗(yàn)證上述方法的有效性，下面進(jìn)行兩組實(shí)驗(yàn)。測(cè)區(qū)位于石家莊市某地，其背景磁化方向I=55°，磁偏角D=-16°，實(shí)測(cè)區(qū)域如圖6所示。實(shí)驗(yàn)所用探頭是利用英國(guó)Bartington公司三軸磁通門(mén)傳感器搭建的十字型磁梯度張量探頭，測(cè)試系統(tǒng)主要包括十字型探頭和數(shù)字采集模塊及軟件操作終端。實(shí)驗(yàn)過(guò)程中探頭固定在無(wú)磁實(shí)驗(yàn)臺(tái)架上，利用掃描方式對(duì)待測(cè)區(qū)域的每一個(gè)觀測(cè)點(diǎn)進(jìn)行測(cè)量。

4.1 孤立磁源

第1組實(shí)驗(yàn)是對(duì)一塊長(zhǎng)方體鐵板進(jìn)行測(cè)試。在地理坐標(biāo)系下，鐵板中心坐標(biāo)為(1.25 m,0.85 m,0.45 m)，尺寸為0.7 m×0.5 m×0.1 m. 獲得的實(shí)測(cè)Azz數(shù)據(jù)如圖7所示。利用NSSRP互相關(guān)估計(jì)法對(duì)其磁化方向進(jìn)行估計(jì)，獲得的估計(jì)結(jié)果如表3所示。從表3中可以看出，獲得的磁化方向估計(jì)值與背景磁場(chǎng)存在較大差異，證明了磁性體具有強(qiáng)剩磁存在。利用NSSRP互相關(guān)法得到的磁化方向估計(jì)值作為總磁化方向，分別利用Azz分量和磁總場(chǎng)模量進(jìn)行反演，獲得的三維反演結(jié)果為磁化強(qiáng)度大于0.3×maxm(α)所顯示的結(jié)果，如圖8和圖9所示。根據(jù)獲得的反演結(jié)果可以看出，利用Azz數(shù)據(jù)進(jìn)行反演能夠獲得更好的橫向和縱向分布結(jié)果，其中在Azz反演過(guò)程中存在兩塊磁性體的原因是反演過(guò)程中受測(cè)量誤差和環(huán)境因素等干擾所導(dǎo)致的。根據(jù)表4中LRSME的計(jì)算結(jié)果也可以看出，利用Azz反演獲得的計(jì)算結(jié)果具有更高精度。

表3 磁化方向估計(jì)結(jié)果

方法LRSMEAzz磁總場(chǎng)模量L2133.861 63 915.2L-curve0.849 71.475 2UPRE0.678 50.838 1

4.2 多目標(biāo)磁源

第2組實(shí)驗(yàn)對(duì)水平圓柱體鐵筒和長(zhǎng)方體鐵塊的組合磁異常進(jìn)行測(cè)試。長(zhǎng)方體鐵塊位于測(cè)區(qū)東南方向，中心坐標(biāo)為(1.5 m，0.4 m，0.3 m)，尺寸為0.23 m×0.36 m×0.20 m，鐵筒位于測(cè)區(qū)西北方向，中心坐標(biāo)為(0.6 m，1.2 m，0.2 m)，圓筒直徑為0.1 m，長(zhǎng)為0.46 m. 獲得的實(shí)測(cè)數(shù)據(jù)如圖10所示。針對(duì)組合鐵磁性體，分別利用Azz分量和磁總場(chǎng)模量反演方法進(jìn)行反演，獲得的三維反演結(jié)果如圖11和圖12所示。

根據(jù)獲得的反演結(jié)果可以看出，利用Azz數(shù)據(jù)進(jìn)行反演時(shí)，利用背景磁場(chǎng)的磁化方向進(jìn)行反演，獲得的反演結(jié)果與磁性體的實(shí)際位置偏差較大，在剩余磁化影響下，對(duì)磁性體的實(shí)際反演結(jié)果產(chǎn)生了較大的影響。利用磁總場(chǎng)模量數(shù)據(jù)進(jìn)行反演時(shí)，能夠獲得更好的橫向和縱向分布結(jié)果，根據(jù)表5中LRSME的計(jì)算結(jié)果也可以看出，該方法獲得的計(jì)算結(jié)果具有更高精度。單次迭代時(shí)，利用Azz進(jìn)行反演的時(shí)間為0.106 s，利用磁總場(chǎng)模量異常進(jìn)行反演的時(shí)間為0.342 s，這是因?yàn)樵谟?jì)算過(guò)程中后者需要3個(gè)靈敏度矩陣，其計(jì)算成本約為前者方法的3倍。

5 結(jié)論

本文提出了強(qiáng)剩磁條件下磁性目標(biāo)三維正則化聚焦反演方法，向目標(biāo)函數(shù)中引入深度加權(quán)矩陣和最小支撐矩陣并利用SVD法對(duì)目標(biāo)函數(shù)進(jìn)行分解，通過(guò)UPRE對(duì)正則化參數(shù)進(jìn)行自適應(yīng)選擇。在反演數(shù)據(jù)的選擇上，針對(duì)孤立磁源，通過(guò)磁梯度張量數(shù)據(jù)進(jìn)行反演；針對(duì)多目標(biāo)磁源，采用磁總場(chǎng)模量數(shù)據(jù)進(jìn)行反演。得到如下主要結(jié)論：

1)利用深度加權(quán)矩陣和最小支撐矩陣對(duì)經(jīng)典Tikhonov正則化理論框架下的反演模型進(jìn)行約束，提高了反演解的穩(wěn)定性。

表5 平均相對(duì)誤差LRSME

2)采用UPRE對(duì)迭代過(guò)程中的正則化參數(shù)進(jìn)行選擇，實(shí)現(xiàn)了正則化參數(shù)選擇的自適應(yīng)功能。

3)提出了適用于多種磁源情況下的磁性目標(biāo)反演方法，能夠同時(shí)對(duì)多種磁性體進(jìn)行反演，并獲得較高的反演精度。