基于形態反演的小尺度磁性目標三維重建方法

尹 剛， 張 林

(中國空氣動力研究與發展中心 高速空氣動力研究所， 四川 綿陽 621000)

小尺度磁性目標的三維重建是以磁性體在測量面產生的磁異常信號為先驗信息，對磁性目標的三維形狀、位置、磁導率(磁化強度)等參數進行求解，進而實現磁性體三維姿態的準確重建.其在重要軍事目標的精確識別[1]、未爆彈銷毀[2]、航天器內部磁性分布研究[3]等領域均具有十分重要的實際應用價值.

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

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

1 模型初始化及最優增長模塊的選擇

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

圖1 形態反演算法流程

在磁性目標三維模型的迭代更新中，需要依據一定的選擇機制選擇最優的增長模塊，為此，建立模塊選擇函數，每次均選擇待增長模塊集中使得模塊選擇函數取值最小的長方體作為最優增長模塊：

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

(1)

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

(2)

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

(3)

由形態反演算法流程可得模型增長示意，如圖2所示，并且在尋找待反演空間中最優解的過程中存在以下約束準則：

(1) 模型增長是連續的，新得到的待增長模塊至少有一個面與當前模型共面，且每一個當前模型都存在一個待增長模塊集，最優增長模塊必須從待增長模塊集中選擇.

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

(3) 整個待反演空間中每個網格的磁化強度只有0和預設值m兩種選擇.

磁性目標形態反演的過程即為反演模型所產生的磁梯度張量場一步步地逼近已知測量數據的過程，并且模型的每次增長都是當前模型最優的，因此，可以認為形態反演所得到的目標模型具有較高的可信度.

圖2 形態反演中模型增長二維示意圖

2 迭代終止準則的建立

在整個形態反演過程中，每次均選擇待增長模塊集中使模塊選擇函數取值最小的長方體作為最優增長模塊以更新反演模型，更新過程是持續的.因此，需要建立適當的目標函數，當目標函數達到事先設定的閾值時，則終止模型的迭代更新，得到的當前模型即為形態反演得到的磁性目標的三維形狀.迭代終止準則的建立方式有多種，其中一種為建立迭代過程中的目標函數，表達式為

(4)

式中：φnew(m)和φnew(m)為加入最優增長模塊后的計算值；φold(m)和φold(m)為加入最優增長模塊前的計算值；τ為比例因子，以調整式(4)前后兩項的重要度.

當目標函數小于事先設定的閾值δ時，則表明增加一個長方體對磁性目標模型產生的磁異常場與測量值之間的擬合度的貢獻較小，此時終止模型的迭代更新.

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

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

圖3 磁性目標在測量平面內產生的磁梯度張量場及張量不變量

3 模擬分析

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

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

圖4 計算得到的傾斜角

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

基于Helbig方法估計得到的磁性目標的磁傾角和磁偏角[7]及不同窗口時計算得到的多個磁傾角和磁偏角方差的倒數如圖5和圖6所示，計算時滑動窗口大小選擇為：3×3～7×7.由圖5和圖6可知，Helbig方法不僅準確估計得到了磁性目標的中心位置，也估計出了磁性目標的磁傾角為20° 和磁偏角為35°.

圖5 磁傾角方差的倒數及磁傾角的平均值

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

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

圖6 磁偏角方差的倒數及磁偏角的平均值

圖7 局部波數法反演得到的磁性目標的垂向位置的分布情況

為清晰地給出形態反演中模型增長的過程，圖8展示了模型增長的前3步、第52～54步以及最后3步，圖中藍色直線表示真實磁性目標的外部框架，黑色箭頭指示出了該迭代步數相對于前一迭代步數時的長方體增長模塊.

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

模型增長到不同階段時對應測量面得到的張量各分量如圖9～12所示.對比不同階段反演得到的張量分量與圖3所示的長方體在測量面產生的理論張量分量可知，形態反演的模型增長過程也是磁梯度張量反演估計值接近真實值的過程.

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

圖8 形態反演中模型增長過程示意圖

圖9 初始模型在測量面產生的磁梯度張量分量

圖10 模型增長到第52步時在測量面產生的磁梯度張量分量

圖11 模型增長到第90步時在測量面產生的磁梯度張量分量

圖12 模型增長到第157步時在測量面產生的磁梯度張量分量

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

圖13 測量面上張量各分量理論值與反演值之間的RMSE隨模型增長步數的變化趨勢

圖14 長方體磁性目標的三維重建結果

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

4 結語

本文針對小尺度磁性目標的三維重建方法開展研究，將地質勘探領域的反演方法擴展到小尺度磁性目標中，建立了適用于小尺度磁性目標的三維重建方法.相比地質勘探領域的磁性參數反演，本文提出的小尺度磁性目標的三維重建方法增加了反演預處理環節，即利用磁梯度張量異常場數據進行待反演空間范圍和目標磁化方向等參數的估計，將反演空間縮小在磁性目標分布范圍內，有效降低了反演多解性、提高了反演精度和反演速度.實驗結果表明：該重建結果與真實目標相似度較高，具有較高的重建可信度和一定的工程應用價值.當然，本文僅考慮了磁性目標均勻磁化的情況，后續將重點研究非均勻磁化的磁性目標體的三維重建問題，并且將構建實際的磁梯度張量測量系統進一步驗證所提重建方法的工程實用性.