尹 剛， 張 林

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 高速空氣動(dòng)力研究所， 四川 綿陽 621000)

小尺度磁性目標(biāo)的三維重建是以磁性體在測(cè)量面產(chǎn)生的磁異常信號(hào)為先驗(yàn)信息，對(duì)磁性目標(biāo)的三維形狀、位置、磁導(dǎo)率(磁化強(qiáng)度)等參數(shù)進(jìn)行求解，進(jìn)而實(shí)現(xiàn)磁性體三維姿態(tài)的準(zhǔn)確重建.其在重要軍事目標(biāo)的精確識(shí)別[1]、未爆彈銷毀[2]、航天器內(nèi)部磁性分布研究[3]等領(lǐng)域均具有十分重要的實(shí)際應(yīng)用價(jià)值.

目前國(guó)內(nèi)外對(duì)小尺度磁性目標(biāo)三維重建技術(shù)的相關(guān)研究仍處于起步階段，類似的磁源反演技術(shù)主要應(yīng)用于地質(zhì)勘探方面，可分為三維形態(tài)反演和物性反演[4].形態(tài)反演是在測(cè)量面下半空間場(chǎng)源體給定磁性參數(shù)大小的基礎(chǔ)上，利用觀測(cè)異常數(shù)據(jù)來擬合幾何體(如多邊形或多面體)形狀，通過幾何體的形態(tài)大小模擬目標(biāo)體的三維姿態(tài).物性反演是將觀測(cè)區(qū)域?qū)?yīng)的測(cè)量面下半空間離散成規(guī)則的長(zhǎng)方體單元，通過反演的方法估計(jì)每個(gè)離散單元的磁性值，由磁性的分布勾繪出場(chǎng)源的三維姿態(tài).形態(tài)和物性反演均可實(shí)現(xiàn)目標(biāo)的三維重建，但物性反演需要求解關(guān)于未知磁性參數(shù)的欠定方程組，存在解不唯一且不穩(wěn)定的問題[5].

盡管小尺度磁性目標(biāo)的三維重建問題可以借鑒地質(zhì)勘探中的磁源反演技術(shù)，但是兩者又存在一定的差別：首先磁性目標(biāo)的三維重建問題對(duì)磁性參數(shù)反演的精度要求更高；其次在參數(shù)反演的基礎(chǔ)上需要勾繪出磁性目標(biāo)的三維姿態(tài).因此，本文在借鑒地質(zhì)勘探領(lǐng)域磁源形態(tài)反演技術(shù)的基礎(chǔ)上，建立適合小尺度磁性目標(biāo)的三維重建技術(shù).利用比磁總場(chǎng)數(shù)據(jù)信息更為豐富的磁梯度張量數(shù)據(jù)[6]進(jìn)行磁性目標(biāo)的三維重建，首先對(duì)測(cè)得的磁梯度張量數(shù)據(jù)進(jìn)行傾斜角、局部波數(shù)和Helbig估計(jì)計(jì)算，得到磁性目標(biāo)水平方向和垂直方向的大致分布范圍和近似磁化方向；然后對(duì)磁性目標(biāo)待反演空間進(jìn)行網(wǎng)格劃分，確定初始模型并建立待增長(zhǎng)模塊集，進(jìn)而在待增長(zhǎng)模塊集中選擇最優(yōu)的增長(zhǎng)模塊對(duì)當(dāng)前模型進(jìn)行迭代更新；通過模型的不斷增長(zhǎng)最終實(shí)現(xiàn)磁性目標(biāo)的三維重建.

1 模型初始化及最優(yōu)增長(zhǎng)模塊的選擇

為得到較為準(zhǔn)確的形態(tài)反演結(jié)果，本文采用不求解反演方程，而是在待反演空間的所有解中找到一個(gè)最優(yōu)解的方法.其基本思想是：首先利用測(cè)得的磁梯度張量數(shù)據(jù)建立磁局部異常特征和磁性目標(biāo)邊界的相互關(guān)系，利用傾斜角[7]和局部波數(shù)[8]計(jì)算，估計(jì)得到磁性目標(biāo)水平和垂向的大致分布范圍，進(jìn)而確定較為緊湊的磁性目標(biāo)的待反演空間范圍并進(jìn)行網(wǎng)格劃分；然后利用Helbig方法[7]估計(jì)磁性目標(biāo)的近似磁化方向，在網(wǎng)格中選定一個(gè)長(zhǎng)方體作為初始模型并對(duì)其磁化強(qiáng)度賦值，建立待增長(zhǎng)模塊集，并使該模塊集中的每一個(gè)長(zhǎng)方體至少有一個(gè)面與當(dāng)前模型中的面重合，進(jìn)而通過一定的選擇機(jī)制在待增長(zhǎng)模塊集中選擇最優(yōu)的增長(zhǎng)模塊對(duì)當(dāng)前模型進(jìn)行迭代更新；通過模型的不斷增長(zhǎng)最終實(shí)現(xiàn)磁性目標(biāo)三維形狀的構(gòu)建.故小尺度磁性目標(biāo)三維重建的算法流程如圖1所示.

圖1 形態(tài)反演算法流程

在磁性目標(biāo)三維模型的迭代更新中，需要依據(jù)一定的選擇機(jī)制選擇最優(yōu)的增長(zhǎng)模塊，為此，建立模塊選擇函數(shù)，每次均選擇待增長(zhǎng)模塊集中使得模塊選擇函數(shù)取值最小的長(zhǎng)方體作為最優(yōu)增長(zhǎng)模塊：

Γ(m)=φ(m)+ρφ(m)+κθ(m)=

(1)

式中：m為長(zhǎng)方體的磁化矢量；N為形態(tài)反演時(shí)使用的磁梯度張量分量的個(gè)數(shù)；ρ,κ為正則因子，用以調(diào)整目標(biāo)函數(shù)中3項(xiàng)的相互比例；φ(m)用相關(guān)性計(jì)算確保反演模型產(chǎn)生的磁梯度張量數(shù)據(jù)與測(cè)得的磁異常數(shù)據(jù)在整個(gè)測(cè)量區(qū)域上的整體相似性，并克服φ(m)約束可能出現(xiàn)的異常點(diǎn)問題，φ(m)為反演模型產(chǎn)生的磁梯度數(shù)據(jù)與測(cè)得的磁梯度數(shù)據(jù)之間的差值，可克服僅僅用φ(m)約束可能存在的整體偏大或整體偏小問題，θ(m)通過計(jì)算長(zhǎng)方體單元之間的距離來約束反演模型，使其具有較好的整體性且不朝某一個(gè)方向無限延伸.φk(m)表示第k個(gè)磁梯度張量分量在整個(gè)測(cè)量區(qū)域上的整體相似性，φk(m)為反演模型產(chǎn)生的第k個(gè)磁梯度張量分量數(shù)據(jù)與測(cè)得的該磁梯度張量分量數(shù)據(jù)之間的差值.具體計(jì)算公式為

(2)

式中：cov(·)表示求協(xié)方差；D(·)表示求方差；dk為在測(cè)量平面內(nèi)測(cè)得的磁梯度張量的第k個(gè)分量的數(shù)據(jù)向量，F(xiàn)k為對(duì)應(yīng)第k個(gè)分量的核函數(shù)矩陣，由單個(gè)長(zhǎng)方體的磁梯度張量正演公式計(jì)算得到[9]；Δx, Δy, Δz為當(dāng)前模型在3個(gè)方向上的寬度；M為當(dāng)前模型所包含的長(zhǎng)方體個(gè)數(shù)；lj為待選擇增長(zhǎng)模塊中心與當(dāng)前模型中第j個(gè)長(zhǎng)方體中心之間的距離；ζk為尺度因子，可由磁梯度張量數(shù)據(jù)的測(cè)量值和正演值求解得到，求解公式為

(3)

由形態(tài)反演算法流程可得模型增長(zhǎng)示意，如圖2所示，并且在尋找待反演空間中最優(yōu)解的過程中存在以下約束準(zhǔn)則：

(1) 模型增長(zhǎng)是連續(xù)的，新得到的待增長(zhǎng)模塊至少有一個(gè)面與當(dāng)前模型共面，且每一個(gè)當(dāng)前模型都存在一個(gè)待增長(zhǎng)模塊集，最優(yōu)增長(zhǎng)模塊必須從待增長(zhǎng)模塊集中選擇.

(2) 待增長(zhǎng)模塊集中包含的每個(gè)模塊均至少有一個(gè)面與當(dāng)前模型共面(如圖2所示)，當(dāng)前模型迭代更新時(shí)，待增長(zhǎng)模塊集也要同時(shí)更新.

(3) 整個(gè)待反演空間中每個(gè)網(wǎng)格的磁化強(qiáng)度只有0和預(yù)設(shè)值m兩種選擇.

磁性目標(biāo)形態(tài)反演的過程即為反演模型所產(chǎn)生的磁梯度張量場(chǎng)一步步地逼近已知測(cè)量數(shù)據(jù)的過程，并且模型的每次增長(zhǎng)都是當(dāng)前模型最優(yōu)的，因此，可以認(rèn)為形態(tài)反演所得到的目標(biāo)模型具有較高的可信度.

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

2 迭代終止準(zhǔn)則的建立

在整個(gè)形態(tài)反演過程中，每次均選擇待增長(zhǎng)模塊集中使模塊選擇函數(shù)取值最小的長(zhǎng)方體作為最優(yōu)增長(zhǎng)模塊以更新反演模型，更新過程是持續(xù)的.因此，需要建立適當(dāng)?shù)哪繕?biāo)函數(shù)，當(dāng)目標(biāo)函數(shù)達(dá)到事先設(shè)定的閾值時(shí)，則終止模型的迭代更新，得到的當(dāng)前模型即為形態(tài)反演得到的磁性目標(biāo)的三維形狀.迭代終止準(zhǔn)則的建立方式有多種，其中一種為建立迭代過程中的目標(biāo)函數(shù)，表達(dá)式為

(4)

式中：φnew(m)和φnew(m)為加入最優(yōu)增長(zhǎng)模塊后的計(jì)算值；φold(m)和φold(m)為加入最優(yōu)增長(zhǎng)模塊前的計(jì)算值；τ為比例因子，以調(diào)整式(4)前后兩項(xiàng)的重要度.

當(dāng)目標(biāo)函數(shù)小于事先設(shè)定的閾值δ時(shí)，則表明增加一個(gè)長(zhǎng)方體對(duì)磁性目標(biāo)模型產(chǎn)生的磁異常場(chǎng)與測(cè)量值之間的擬合度的貢獻(xiàn)較小，此時(shí)終止模型的迭代更新.

迭代終止的判斷還可以利用測(cè)量面上磁梯度張量各分量真實(shí)值與反演估計(jì)值之間的方均根誤差(RMSE)隨模型增長(zhǎng)步數(shù)的變化規(guī)律來實(shí)現(xiàn)[10].當(dāng)反演過程從初始模型向真實(shí)物體接近時(shí)，隨著模型增長(zhǎng)步數(shù)的增加，各張量分量的RMSE逐漸變小，直至反演模型最接近真實(shí)物體時(shí)，各張量分量的RMSE達(dá)到最小值.此時(shí)若繼續(xù)進(jìn)行模型增長(zhǎng)計(jì)算，各分量的RMSE將逐步增大.因此，在模型增長(zhǎng)過程中，磁梯度張量分量真實(shí)值與反演估計(jì)值之間的RMSE極小值點(diǎn)可作為模型迭代增長(zhǎng)的終止點(diǎn)，其對(duì)應(yīng)的反演模型即為形態(tài)反演得到的磁性目標(biāo)的三維形狀.

另外，當(dāng)實(shí)際測(cè)量數(shù)據(jù)中存在較大的測(cè)量噪聲時(shí)，將導(dǎo)致ξ(m)的取值始終都較大且張量分量估計(jì)值與理論值之間的RMSE隨步數(shù)增大的變化規(guī)律不明顯.此時(shí)，可結(jié)合傾斜角計(jì)算得到的水平邊界約束范圍進(jìn)行判斷，若估計(jì)得到的磁性目標(biāo)模型已遠(yuǎn)遠(yuǎn)超出傾斜角估計(jì)得到的水平邊界范圍，則認(rèn)為測(cè)量數(shù)據(jù)中噪聲較大導(dǎo)致閾值δ較小，需終止迭代過程重新選擇δ進(jìn)行形態(tài)反演.

圖3 磁性目標(biāo)在測(cè)量平面內(nèi)產(chǎn)生的磁梯度張量場(chǎng)及張量不變量

3 模擬分析

為驗(yàn)證所提形態(tài)反演方法用于磁性目標(biāo)三維重建中的有效性，以空間中存在的磁性長(zhǎng)方體為例進(jìn)行數(shù)值模擬實(shí)驗(yàn).其中，長(zhǎng)方體中心坐標(biāo)為(0,0,10 m)，邊長(zhǎng)分別為12，12和8 m，考慮較為一般性的情況，設(shè)置磁化強(qiáng)度為40 A/m，磁化方向的傾角和偏角分別為20° 和35°.建立z方向豎直向下為正的右手坐標(biāo)系，考慮實(shí)際中運(yùn)動(dòng)磁梯度張量系統(tǒng)的測(cè)量條件，設(shè)置水平x和y方向上的采樣間隔均為2 m，則在z=-1 m平面測(cè)得的磁性目標(biāo)產(chǎn)生的磁梯度張量的6個(gè)分量Bxx、Bxy、Bxz、Byy、Byz、Bzz和計(jì)算得到的張量不變量(Normalized Source Strength, NSS)[11]如圖3所示.

計(jì)算得到的傾斜角[7]如圖4所示，由于該磁性目標(biāo)的磁傾角較小，計(jì)算的得到的傾斜角I1和I2均出現(xiàn)正負(fù)伴生的現(xiàn)象，而I3并未出現(xiàn)此種情況，所以，模擬實(shí)驗(yàn)以傾斜角I3的識(shí)別結(jié)果進(jìn)行目標(biāo)反演空間的水平約束，取水平方向上x和y的約束范圍均為 -8～8 m.

圖4 計(jì)算得到的傾斜角

值得注意的是，在此得到的磁性目標(biāo)反演空間的水平約束并非磁性目標(biāo)的水平邊界，而是要大于目標(biāo)的水平邊界，一方面是為了避免水平邊界識(shí)別存在的誤差導(dǎo)致反演空間較小，另一方面是為了給待反演空間劃分長(zhǎng)方體單元模型預(yù)留空間余量.另外，盡管磁梯度張量異常數(shù)據(jù)的水平方向測(cè)量范圍為 -20～20 m，但通過傾斜角識(shí)別得到磁性目標(biāo)存在的近似水平區(qū)域，將反演空間范圍大幅度減小，在提高反演效率的基礎(chǔ)上也可以有效提高反演精度.

基于Helbig方法估計(jì)得到的磁性目標(biāo)的磁傾角和磁偏角[7]及不同窗口時(shí)計(jì)算得到的多個(gè)磁傾角和磁偏角方差的倒數(shù)如圖5和圖6所示，計(jì)算時(shí)滑動(dòng)窗口大小選擇為：3×3～7×7.由圖5和圖6可知，Helbig方法不僅準(zhǔn)確估計(jì)得到了磁性目標(biāo)的中心位置，也估計(jì)出了磁性目標(biāo)的磁傾角為20° 和磁偏角為35°.

圖5 磁傾角方差的倒數(shù)及磁傾角的平均值

利用局部波數(shù)[8]計(jì)算得到的磁性目標(biāo)的垂向位置的分布情況如圖7所示，圖中直方圖顏色僅用于區(qū)分不同深度的估計(jì)個(gè)數(shù)，計(jì)算時(shí)滑動(dòng)窗口大小選擇為6×6.由圖7可知，盡管反演得到的磁性目標(biāo)垂向分布位置跨度較大，約為8～20 m，但其在每個(gè)深度的分布情況各有不同，從左圖中可看出該目標(biāo)大約集中分布在8～14 m，由右圖中歸一化直方圖可得，在深度為14 m處時(shí)直方圖總值可達(dá)到90%，為此，本實(shí)驗(yàn)選擇垂向約束范圍為0～20 m，磁性目標(biāo)垂向中心位置為14 m.

選擇每個(gè)長(zhǎng)方體的大小均為2 m×2 m×2 m，然后結(jié)合圖4和圖7得到的反演空間范圍約束及圖5和圖6得到的磁化方向估計(jì)結(jié)果，設(shè)定形態(tài)反演時(shí)初始長(zhǎng)方體模塊的中心位置為(-1 m，-1 m，9 m)，磁化方向中的磁傾角為20°，磁偏角為35°，進(jìn)而利用式(1)所示的最優(yōu)模塊增長(zhǎng)函數(shù)及式(4)所示的迭代終止準(zhǔn)則進(jìn)行長(zhǎng)方體的形態(tài)反演.

圖6 磁偏角方差的倒數(shù)及磁偏角的平均值

圖7 局部波數(shù)法反演得到的磁性目標(biāo)的垂向位置的分布情況

為清晰地給出形態(tài)反演中模型增長(zhǎng)的過程，圖8展示了模型增長(zhǎng)的前3步、第52～54步以及最后3步，圖中藍(lán)色直線表示真實(shí)磁性目標(biāo)的外部框架，黑色箭頭指示出了該迭代步數(shù)相對(duì)于前一迭代步數(shù)時(shí)的長(zhǎng)方體增長(zhǎng)模塊.

由圖可知，當(dāng)初始模型給定后，其對(duì)應(yīng)的待增長(zhǎng)模塊集也是確定的，后續(xù)的每次迭代，算法均自動(dòng)從待增長(zhǎng)模塊集中選擇最優(yōu)的增長(zhǎng)模塊添加到現(xiàn)有模型中，同時(shí)更新待增長(zhǎng)模塊集為下一次迭代做準(zhǔn)備.整個(gè)迭代過程中，最優(yōu)增長(zhǎng)模塊的選擇一直遵循式(1)所示的準(zhǔn)則，在待增長(zhǎng)模塊集中尋找函數(shù)最小值對(duì)應(yīng)的模塊單元，使得每次都把待增長(zhǎng)模塊集中對(duì)反演貢獻(xiàn)最大的長(zhǎng)方體模塊加入到當(dāng)前模型中.因此，整個(gè)迭代過程即為一個(gè)優(yōu)化過程，通過每一步都選擇最優(yōu)的增長(zhǎng)模塊進(jìn)而實(shí)現(xiàn)最優(yōu)的磁性目標(biāo)形態(tài)反演.

模型增長(zhǎng)到不同階段時(shí)對(duì)應(yīng)測(cè)量面得到的張量各分量如圖9～12所示.對(duì)比不同階段反演得到的張量分量與圖3所示的長(zhǎng)方體在測(cè)量面產(chǎn)生的理論張量分量可知，形態(tài)反演的模型增長(zhǎng)過程也是磁梯度張量反演估計(jì)值接近真實(shí)值的過程.

分析模擬實(shí)驗(yàn)過程可知，由于式(1)中φ(m)的約束，在模型增長(zhǎng)中的每一步，測(cè)量面得到的磁梯度張量分量反演估計(jì)值與理論值均有較高的分布相似度；由于φ(m)的約束，模型朝著反演估計(jì)值與理論值數(shù)值差異變小的方向增長(zhǎng)；由于θ(m)的約束，保證了反演模型不朝垂直向下方向無限延伸，所以，形態(tài)反演通過模型的不斷增長(zhǎng)實(shí)現(xiàn)了磁性目標(biāo)三維形狀的重建.

圖8 形態(tài)反演中模型增長(zhǎng)過程示意圖

圖9 初始模型在測(cè)量面產(chǎn)生的磁梯度張量分量

圖10 模型增長(zhǎng)到第52步時(shí)在測(cè)量面產(chǎn)生的磁梯度張量分量

圖11 模型增長(zhǎng)到第90步時(shí)在測(cè)量面產(chǎn)生的磁梯度張量分量

圖12 模型增長(zhǎng)到第157步時(shí)在測(cè)量面產(chǎn)生的磁梯度張量分量

隨著模型增長(zhǎng)步數(shù)的增加，測(cè)量面上磁梯度張量各分量真實(shí)值與反演估計(jì)值之間的RMSE變化趨勢(shì)如圖13所示，該圖也驗(yàn)證了形態(tài)反演的過程是張量分量反演估計(jì)值向真實(shí)值接近的過程這一初始算法思路，其中RMSE在模型增長(zhǎng)到第157步時(shí)達(dá)到極小值，并且此時(shí)張量各分量估計(jì)值與真實(shí)值的RMSE接近0.因此，可選擇該參數(shù)作為迭代增長(zhǎng)的終止步數(shù)，將此增長(zhǎng)步數(shù)下得到的反演模型作為磁性目標(biāo)形態(tài)反演的最終結(jié)果.將其真實(shí)長(zhǎng)方體的三維形狀對(duì)比可得：所提形態(tài)反演方法得到的三維重建結(jié)果具有較高的可信度.當(dāng)然，由于磁梯度張量數(shù)據(jù)的快速衰減特性，重建得到的磁性目標(biāo)的縱向分辨率略差于橫向分辨率.

圖13 測(cè)量面上張量各分量理論值與反演值之間的RMSE隨模型增長(zhǎng)步數(shù)的變化趨勢(shì)

圖14 長(zhǎng)方體磁性目標(biāo)的三維重建結(jié)果

為了得到直觀的可視效果并反映小尺度磁性目標(biāo)的三維特征，利用上述的反演結(jié)果進(jìn)行三維模型的構(gòu)建，三維重建結(jié)果如圖14所示，盡管重建結(jié)果在長(zhǎng)方體的棱角處存在一定的偏差，但重建結(jié)果與真實(shí)磁性目標(biāo)整體上具有較高的相似度，可用于小尺度磁性目標(biāo)的識(shí)別.后續(xù)可通過更為精細(xì)的網(wǎng)格劃分，選擇更小的模型增長(zhǎng)模塊，長(zhǎng)方體棱角處的反演將更為精確.另外，值得說明的是，本文所提方法不僅僅局限于長(zhǎng)方體狀磁性目標(biāo)的三維重建，可適用于任何形狀的磁性目標(biāo)體.

4 結(jié)語

本文針對(duì)小尺度磁性目標(biāo)的三維重建方法開展研究，將地質(zhì)勘探領(lǐng)域的反演方法擴(kuò)展到小尺度磁性目標(biāo)中，建立了適用于小尺度磁性目標(biāo)的三維重建方法.相比地質(zhì)勘探領(lǐng)域的磁性參數(shù)反演，本文提出的小尺度磁性目標(biāo)的三維重建方法增加了反演預(yù)處理環(huán)節(jié)，即利用磁梯度張量異常場(chǎng)數(shù)據(jù)進(jìn)行待反演空間范圍和目標(biāo)磁化方向等參數(shù)的估計(jì)，將反演空間縮小在磁性目標(biāo)分布范圍內(nèi)，有效降低了反演多解性、提高了反演精度和反演速度.實(shí)驗(yàn)結(jié)果表明：該重建結(jié)果與真實(shí)目標(biāo)相似度較高，具有較高的重建可信度和一定的工程應(yīng)用價(jià)值.當(dāng)然，本文僅考慮了磁性目標(biāo)均勻磁化的情況，后續(xù)將重點(diǎn)研究非均勻磁化的磁性目標(biāo)體的三維重建問題，并且將構(gòu)建實(shí)際的磁梯度張量測(cè)量系統(tǒng)進(jìn)一步驗(yàn)證所提重建方法的工程實(shí)用性.