鄭建擁 范紅波 張 琪 李志寧
(①陸軍工程大學石家莊校區,河北石家莊 050003; ②中國人民解放軍94019部隊,新疆和田 848000 )
利用地磁場異常探測地下或水下物體,一直是勘探地球物理界的研究熱點。目前地磁技術研究的主要對象是地下或水下大型異常體(如潛艇),而對諸如管道、未爆彈藥等小型鐵磁性物體,始終缺乏成熟的技術進行探測。當前,磁梯度張量方法受到廣泛關注,成為探測小型磁性目標體形狀、大小、姿態等信息的新手段。
孫刃[1]較系統地分析總結了利用磁異常正反演識別目標的方法,并將其用于水下沉船、管線、未爆軍火等目標的探測;吳國超[2]提出利用張量矩陣特征值和6個矩陣元素對大型地質體進行邊界識別,分辨率較高,試驗效果很好;王林飛等[3]利用磁梯度張量研究大型磁源體的位置與邊界;李金朋等[4]研究了地下小目標體的反演及三維成像;尹剛等[5]利用張量不變量和矩陣特征值定位磁目標體;朱慧慧等[6]建立了基于磁偶極子構造法幾何建模的地下管線磁異常正演模型,分析了多種磁化因素對鐵管磁異常正演結果的影響,為高精度地下管線磁探測奠定了理論基礎。上述方法主要針對磁性目標的二維邊界識別、三維姿態反演、磁性目標定位等,是局部磁異常識別方法,而且無法避免磁測數據的正、反演,運算量大,受測量數據精度的限制較大,因而不能清晰地表現磁測目標體。謝永茂[7]提出了建立樣本庫,通過模板匹配識別地下磁性體,提供了一種新的磁異常體識別思路。
本文首先建立地下小型磁性目標體不同形狀、不同方位條件下模型的正演數據庫,分析了9個磁梯度張量矩陣屬性參量,論證其作為分類屬性的可行性; 利用這9個矩陣屬性參量融合識別磁源目標體的形狀,可在一定程度上彌補磁測數據精度的不足; 再將這9個屬性參量用于支持向量機(Support Vector Machine,SVM)對目標進行分類; 然后利用量子粒子群算法優化的支持向量機對仿真模型數據分類,推斷目標體形狀,取得了很好的識別效果; 最后設計實驗驗證其實用性。
地磁總場包括地磁場和異常場,其中地磁場在空間上的變化很小,在小范圍內可以看做一個穩定場,所以磁梯度張量可以認為是磁異常矢量在x、y、z三個方向上的變化率。梯度張量元素受地磁場傾向、偏角影響較小,通過反演能夠更精確地描述場源體的磁化方向和幾何形態[8]。
地磁場是矢量場,由矢量場三分量推導梯度張量矩陣G表達式為
(1)
式中U為磁標勢,矩陣中的9個元素即為磁場矢量B在空間三個方向分量(Bx、By、Bz)的變化率。磁法勘探中地磁場及鐵磁物質產生的異常場可看作無源的靜磁場[5]。由式(1)可知梯度張量矩陣是對稱的,且旋度和散度均為0,因此梯度張量矩陣G的9個元素中,只有5個(Bxx、Byy、Bxy、Bxz、Byz)相互獨立[9]。
作為一種機器學習方法,支持向量機需對分類對象進行特征提取,通過對特征進行分類來識別目標。每一個特征都與分類對象本身相關,且不同對象的特征之間都應該保持離散。特征的數量越多,分類越準確。為保證信息不冗余,不同的特征之間必須互不相關,或相關度很低,這些特征即可作為支持向量機分類的屬性參量。從樣本庫中提取不同屬性參量的數據構造特征向量,作為支持向量機的訓練和測試樣本數據。
選擇磁總場模量(TMI)、矩陣特征值的最大值(ME)、兩個矩陣不變量(I1、I2)作為分類的屬性參量。
TMI可由磁場三分量計算
(2)
設A為n階矩陣,λ為常數,X為n維非零向量。若AX=λX,則λ稱為矩陣A的特征值。因為式(1)是3階張量矩陣,因此最多只有3個特征值。其中最大特征值對磁目標體邊界和形狀都很敏感,識別能力較強。
矩陣不變量能夠表征一個矩陣的特征,經由梯度張量矩陣計算得到的張量不變量不需要額外的處理就可以很好地描述磁場源。該矩陣的不變量蘊含著很多磁目標體的特征信息,可以作為地下磁目標體識別的分類屬性參量。式(1)中矩陣的特征值即是基本不變量[10],矩陣中的其他不變量可由矩陣特征值表示。
定義λ1、λ2、λ3為張量矩陣G的特征值,則其余兩個矩陣不變量為
(3)
除以上4個屬性參量之外,還提取了式(1)中的Bxx、Bxy、Bxz、Byy、Byz等5個獨立元素作為支持向量機分類的屬性參量,總共使用這9個屬性參量識別目標。
支持向量機可用于模式分類和非線性回歸,具有結構風險最小化的優勢,具有很好的泛化性能。支持向量機的算法公式為
(4)

加入松弛因子ξ(i),用于衡量實際指示值Li與支持向量機輸出之間的距離,支持向量機優化的目標函數為[11]
Li(wTXi+b)≥1,i=1,2,…,l
(5)
式中:w即是對目標進行分類的超平面的法向量;l為分類時特征向量的維數;c是懲罰因子。
影響SVM的因素非常多,如懲罰因子c以及核函數參數σ,對于小型樣本分類精準度的影響最大[12]。因此,在應用SVM做分類預測時,這兩個參數的選取尤為關鍵。
一般通過網格劃分來尋找參數c和σ的最佳值。雖然采用網格搜索能夠獲得最高的分類準確率,即全局最優解,但如果要在更大范圍內尋找參數c和σ的最佳值會很費時。本文應用一種通過尋找最優解來求解最佳參數的算法,即量子粒子群優化算法(Quantum Particle Swarm Optimization,QPSO)[11]。
粒子群優化算法是智能計算領域中,除蟻群算法之外的另一種基于群體智能的優化算法。該方法通過粒子在解空間的搜索來尋找最優解。相比基因等其他算法,粒子群優化算法無需交叉、變異、選擇等流程,正確率比較高,全局最優值收斂速度快,參數尋優耗時短。
Sun等[13]認為粒子具有量子性,從量子力學的角度提出了一種新的粒子群優化算法模型,并經實驗證明其具有更好的收斂性。粒子群優化算法[14]具體過程如下。
在有M個粒子進行尋優的J維搜索空間中,假設第i個粒子的位置為Xi=(Xi,1,Xi,2,…,Xi,J),這里t為進化代數; 其每代最優位置為Pi=(pi,1,pi,2,…,pi,J),整個量子粒子群目前所搜索到的全局最優解為Pg=(pg,1,pg,2,…,pg,J)。通過200次迭代,選取全局最優的c和σ。
粒子i的位置為
(6)

最后得到進化t代之后的pi,j,即可得到最優參數c和σ。基于量子粒子群算法優化后的c和σ,可建立最優化QPSO-SVM分類模型。
為了進行仿真驗證,首先構建一個垂直磁化的磁性體進行正演。在觀測面上建立以觀測原點為中心的61×61個網格節點的數據平面。數據平面上每個網格節點都包含該位置的梯度張量信息,即一個由式(1)計算的3階梯度張量矩陣。
設置8組長方體磁性體正演模型,模型中心位置均處于坐標原點,埋深均為5m,磁化率相同,姿態都為水平,但他們的形狀和磁化方向不同,即構建8種不同尺寸的長方體模型。每個模型的磁化偏角范圍為0°~180°,間隔設置為10°(0°與180°重合),這樣得到共180種不同磁化方向組合,總共8×9×180個包含61×61個網格點的數據平面。分別計算平面上61×61個點的TMI、3個不變量和5個張量矩陣元素,這樣每個平面都可以衍生出9個61維的獨立屬性矩陣。具體的數據結構參見圖1。
圖1為9個屬性參量的等值線圖,每幅圖都代表一個屬性參量的屬性矩陣。圖中每個像素點的數值就是矩陣中相應位置的值。
為論證9個屬性參量作為支持向量機分類的屬性參量的可行性,以長、寬、高均為10m、埋深為15m的正方體模型為例進行正演。設磁傾角為90°,磁偏角為0°,磁化率為0.4。
由圖1可見,TMI的等值線分布(圖1a)比較均勻,但磁化體邊角特征不明顯;磁性體矩陣的最大特征值ME(圖1b)和I1(圖1c)相較于其他屬性參量,能更詳細和清晰地描述磁性體的邊界和形狀特征;I2(圖1d)對邊角特征更敏感,能凸顯淺部源; 5個矩陣獨立元素(圖1e~圖1d)中,Bxx描述y軸方向的邊界情況,出現了三個極值點,Bxy可描述磁性體的邊角,Bxz可描述y軸方向的邊界情況,出現了兩個極值點,Byy可描述x軸方向的邊界情況,出現了三個極值點,Byz主要描述x軸方向的邊界,出現了兩個極值點[15]。可見這5個矩陣獨立元素也能較清晰地表現磁性體的部分特征。可以看出,這5個屬性參量都與磁性體的形狀有著較高的相關性,可從不同角度反映磁性體的形狀特征。
圖2為所有樣本的9個屬性參量的分維屬性圖,可據此分析這9個獨立屬性參量對不同形狀和傾角的敏感性。圖中每個點線代表不同的磁偏角和磁傾角。對于傾角的敏感性,可根據每條點線的傾斜角度進行分析,因為每個點都代表不同的傾角,而這些點的排列是有規律的,因此對傾角變化敏感性強的屬性參量就是可用于研究磁化體形狀和姿態特征的分類。

圖1 9個屬性參量的等值線圖
不失一般性,可以認為這9個屬性參量與磁源的形狀、磁化方向等信息有著很高的相關性。由圖2可見,不同形狀時正演結果變化明顯,如基于屬性參量TMI(圖2a)與屬性參量I2(圖2d)可以很明顯地區分出8種形態;而屬性參量I1(圖2c)和Bxx(圖2e)對于磁化角度的敏感性隨著磁偏角的變化明顯不同,分別出現遞增和遞減兩種變化趨勢。從圖2還可以看出,分類屬性值與偏角呈明顯的線性關系。屬性參量Bxx與Byy、Bxz與Byz形態雖然相似,但數值大小相差兩個數量級。由式(1)可知,Bxx與Byy、Bxz與Byz相互獨立,因此它們之間具有離散性,即各屬性參量之間的數值相差很大,遠遠大于屬性參量內部的數值差,不同參量的數值在數域上明顯分別開。因此,可以參與分類。
上述分析表明,這9個屬性參量對磁性體有良好的辨識能力,能夠作為支持向量機分類的屬性參量。每個屬性參量的樣本數據都可以構造一個特征向量,作為支持向量機分類的輸入向量。

圖2 9個屬性參量的分維屬性
從樣本庫中調取數據,計算每個獨立屬性矩陣的最大特征值,分別對應9個特征向量中的某一個元素。對于這8種形狀的磁化體,每種磁化體共有180個屬性參量平面,即可構造出8組9×180的特征矩陣。以此作為PSO-SVM的訓練樣本,并隨機抽選每種形狀下各90組數據為訓練樣本,另外90組數據可作為測試樣本。
在支持向量機中, 8種長方體形狀分別對應的輸出代碼為整數1~8,即類別標簽設置為1~8,分類結果通過數字體現。
在仿真數據輸入前,必須對數據進行歸一化預處理
(7)
式中p、q∈R。歸一化的目的是將原始數據規一化到[0,1]范圍。結果輸出后再進行反歸一化,最后經過處理獲得1~8的分類標簽。
設置量子粒子群優化的初始搜索參數為:最大進化數量為200,種群最大數量為20。經過200次進化尋優,可以得到參數c和σ的最優解分別為77.4375和6.2357。
適應度用于評價粒子尋找最優解的能力強弱,適應度越大粒子尋優能力越強,反之適應度越小。適應度度量函數因具體情況而異。適應度的度量公式為
(8)
圖3為最終適應度曲線,其中最佳適應性曲線為每代粒子中的最佳適應度,平均適應度曲線為每代粒子的平均適應度,其最佳適應度最終為87.36%。

圖3 最終適應度曲線
利用得到的c和σ的最優解構建基于QPSO-SVM的支持向量機。輸入720個訓練樣本,對支持向量機進行訓練,最后用測試樣本進行測試,結果見圖4,識別正確率約為93%。

圖4 測試分類與實際分類對比圖

圖5 實驗系統設計
為驗證本文方法的可行性,采用4個磁通門傳感器構建了平面十字形磁梯度張量系統。圖5為實驗所使用的實驗臺,圖6為實測結果。
實驗采用直徑分別為10、20cm、對應長度為100、70cm的兩種鐵圓柱和一個半徑為20cm、厚度為10cm的鐵圓板三種樣品作為識別對象。分類標簽分別設置為1、2、3。設置水平埋深為2m。搭建邊長(x,y軸)為2m的水平正方形實驗測量系統,x和y軸方向每隔0.1m設置一測量點,測量三種樣品在同一空間位置時的正方形實驗面中20×20個測量點的磁場數據,因此得到一個20×20的矩陣。對三種樣品的4種姿態(與x軸分別成0°、45°、90°和135°)各采集了5組磁異常數據,并將樣本數據以3∶2的比例隨機分成訓練樣本和測試樣本,建立磁異常數據庫。實驗測量前,先測量背景場,對儀器進行校正。

圖6 實測磁異常場分布
以長度為1m,直徑為10cm的鐵質圓柱體為例分析結果。將其放置在實驗臺的正中位置,與實驗臺x軸夾角為45°,測量其地磁總場后減去背景場,即可得到異常場。由圖6可見,實測數據包含很強的干擾和噪聲,主要來自大地、實驗儀器和其他地下大型鐵磁性物體,干擾了對目標體方位、輪廓、邊緣等細節的分析。實驗地點處于中緯度地區,磁傾角約為55°,磁偏角約為-9°,斜磁化使磁異常信號畸變和偏移,提高了信號分析的難度。但是根據圖6的磁異常分布仍然可以較清晰地分辨出目標的方位,也可大致辨別異常體的部分輪廓信息,且表現出與圖1理論仿真信號相似的形態特征。
圖6中屬性參量TMI、I1、I2和ME均在一定程度上包含目標體的形狀和姿態信息。對比圖1的仿真結果可以看出,這4個參量的分布形態雖然發生了偏移和扭曲,但是等值線的形態仍與目標體保持了一定的相似性,可反映出目標的部分輪廓。而另外的5個張量矩陣獨立元素受到斜磁化影響較大,不能直觀地反映目標形態,但是等值線峰值的連線卻可大致反映出目標體的姿態。以上9個屬性參量分布圖表明,通過實測磁異常信號并不能直觀地識別目標,但是實際信號與目標之間仍保持一定的相關性。因此,不失一般性地,認為這9個屬性參量可以用于實際目標分類。這就是本方法的基本思想,對于不能直觀反映目標信息的實測磁場,利用支持向量機對信號的特征進行分類,以達到識別目標的目的。
最后采用PSO -SVM進行分類,結果見圖7,識別正確率約為87.5%。
由實驗可知,本方法對小型鐵磁性物體的識別取得了一定的效果。由于地磁背景場比較強,異常體或目標體的磁異常信號易被背景磁場淹沒,降低了正確率。如圖7所示,識別錯誤主要集中在標簽1和標簽2,有兩個標簽1被識別為2,一個標簽2被識別為1,這是因為這兩種識別對象的形狀相似度較高。因此關于細節識別有待進一步的研究。

圖7 實驗分類結果
(1)利用矩陣最大特征值,構造模式識別的特征向量,分析并提取張量矩陣的9個屬性參量聯合識別磁性目標體,一定程度上彌補了磁測數據精度不足的缺點,也避免了大量磁測數據的計算。
(2)將機器學習的方法引入地下磁性目標體的識別。利用量子粒子群改進的支持向量機識別地下小型目標體的形狀,取得了較好的效果。
(3)本文僅研究磁性異常體形狀的識別,但此方法還可以識別磁目標的姿態、深度等。由實驗結果可知,本文方法對磁異常信號的測量和特征分析仍有很大的研究空間。