黃雙龍，邱 景，*，孫赫軒，張順雨，蘭 天

（1.重慶大學(xué) 光電工程學(xué)院，重慶 400044；2.光電技術(shù)及系統(tǒng)教育部重點實驗室，重慶 400044）

0 引言

當(dāng)前水下成像技術(shù)主要是利用光學(xué)或聲學(xué)等手段在水下獲取目標(biāo)或環(huán)境的圖像信息。水下光學(xué)成像可以提供高分辨率、符合人眼視覺特征的圖像，但受限于能見度，特別是在渾濁的水體中，光學(xué)成像的有效距離很短，難以實現(xiàn)大范圍的海底探測。聲波在水中的衰減遠(yuǎn)比電磁波小，可以實現(xiàn)超遠(yuǎn)距離的水下探測和通信。但聲波的波長與頻率成反比，高頻聲波具有較高的分辨率，而衰減較快，低頻聲波具有較遠(yuǎn)的傳播距離，而分辨率較低，難以同時滿足高分辨率和遠(yuǎn)距離的要求。由于水下目標(biāo)（管道、線纜、沉船等）大多數(shù)帶有磁性，而磁異常信號具有不易受海水、泥沙等介質(zhì)影響的優(yōu)點，適合用于掩埋物的探測，因此磁學(xué)成像技術(shù)被越來越多地應(yīng)用于探測水下目標(biāo)[1-3]。

磁學(xué)成像技術(shù)主要分為磁化率成像法和相關(guān)成像法。磁化率成像法可以根據(jù)磁場數(shù)據(jù)來推斷水下磁性體的位置、形狀、大小和磁化強度等大量信息，有利于識別水下目標(biāo)的類型。但是在缺乏先驗信息的條件下，磁化率成像算法計算速度慢，可能陷入局部最優(yōu)解或收斂困難[4-6]。而相關(guān)成像技術(shù)則是計算實測磁異常與地下各點產(chǎn)生的磁異常的歸一化互相關(guān)系數(shù)，并用相關(guān)系數(shù)的極大值位置估計磁性體的中心位置，過程簡單、穩(wěn)定、易于實現(xiàn)[7]。因此相關(guān)成像技術(shù)可以作為磁化率成像的前置步驟，提供磁性體的位置先驗信息，改善磁化率成像效果。

MAURIELLO等首次將相關(guān)成像法用于磁性目標(biāo)的定位，通過應(yīng)用Schwarz不等式，導(dǎo)出了用于成像磁異常源分布的發(fā)生概率函數(shù)[8]。由于完全忽略先驗信息，相關(guān)成像方法的縱向分辨率較差，除了簡單和孤立的磁性體外，它無法區(qū)分分布在復(fù)雜或緊密間隔情況下的磁性體，因此對于不同深度的多磁性目標(biāo)，相關(guān)成像的效果不佳[4，9]。目前已經(jīng)有許多研究對相關(guān)成像進行了改善，GUO等提出了一種新的磁場總異常垂直梯度相關(guān)成像方法，證明了垂直梯度數(shù)據(jù)的三維相關(guān)成像比單獨使用磁場總異常的分辨率更高[10]。馬國慶等改進了相關(guān)成像技術(shù)，減小了磁性體的實際形狀與通常使用的球體形狀差別較大時所產(chǎn)生的計算誤差，但難以對疊加場進行處理[11]。李金鵬等人提出了基于磁梯度張量的地下小目標(biāo)相關(guān)成像方法，通過迭代計算將概率密度分布轉(zhuǎn)化為實際的物性參數(shù)分布[12]。SUN等從相關(guān)成像結(jié)果中提取空間變化，然后將其轉(zhuǎn)換為空間加權(quán)函數(shù)添加到模型目標(biāo)函數(shù)中，獲得了具有更高可靠性的結(jié)果[13]。MA等利用磁異常數(shù)據(jù)的局部波數(shù)和假設(shè)源生成的合成數(shù)據(jù)的局部波數(shù)的相關(guān)系數(shù)來計算磁源的深度[14]。馬國慶等接著提出了重磁不同階梯度比值相關(guān)成像法，該方法可以準(zhǔn)確計算出磁性體的中心位置，且在二階及以上導(dǎo)數(shù)計算時采用Laplace方程來完成，具有良好的抗噪性[15]。林濤等在均值濾波和插值切割分離的位場異常上，使用相關(guān)成像法對穩(wěn)定向下延拓得到不同深度層上重力及梯度異常進行處理，提高了縱、橫向成像分辨率[16]。丁亞洲推導(dǎo)出了不同階次解析信號垂向?qū)?shù)的相關(guān)成像法，也有效提高了相關(guān)成像的分辨率[17]?？傊?，相關(guān)成像技術(shù)作為一種相對新型的快速自動解釋方法，這些年來許多研究都在不斷改善其成像效果、逐漸提高成像分辨率。

為了進一步提升對不同深度的多磁性目標(biāo)的相關(guān)成像效果，本文使用快速傅里葉變換對測量平面的磁梯度數(shù)據(jù)進行處理，根據(jù)不同深度磁性體的磁梯度數(shù)據(jù)頻譜特性進行自適應(yīng)濾波，將濾波后的數(shù)據(jù)用于成像，有助于提升成像結(jié)果的分辨率，進一步發(fā)揮磁學(xué)成像技術(shù)在探測水下目標(biāo)領(lǐng)域的應(yīng)用。

1 磁梯度數(shù)據(jù)三維相關(guān)成像法

磁梯度張量G是磁感應(yīng)強度B的3個分量各自的空間變化率。在三維直角坐標(biāo)系下分別對磁感應(yīng)強度B的3個分量沿X、Y、Z這3個方向求導(dǎo)，就能得到磁梯度張量G，其表達式為

式中：Bx、By，、Bz為磁感應(yīng)強度B在X、Y、Z這3個方向上的分量；Bij（i，j=x，y，z）為磁梯度張量的各個分量。

由于地磁場通常可以看作無源靜磁場，所以磁感應(yīng)強度B的散度以及旋度都為0，表達式為

因此，磁梯度張量的各個分量之間并非相互獨立，而是滿足

根據(jù)式（3）可以得到，磁梯度張量雖然有9個分量，但其中僅有5個是獨立分量，而且其構(gòu)成的磁梯度張量矩陣是一個對稱矩陣。

與傳統(tǒng)的磁總場異常相比，磁梯度張量能提供更加豐富的異常信息，能更好地反映水下磁性體的細(xì)節(jié)特征[18]。此外，使用磁梯度張量的每個分量進行相關(guān)成像的效果也不盡相同，通常Bzz分量在深度方向上具有最好的成像效果[19]。由于相關(guān)成像具有橫向分辨率較高、縱向分辨率差的特點[20]，本文采用Bzz分量進行相關(guān)成像。

磁梯度數(shù)據(jù)三維相關(guān)成像法的基本流程：1）將待成像的空間劃分成M個均勻的三維網(wǎng)格；2）計算每一個網(wǎng)格可能存在的磁性體在觀測平面引起的磁梯度數(shù)據(jù)；3）計算實測數(shù)據(jù)與計算數(shù)據(jù)的相關(guān)系數(shù)；4）將相關(guān)系數(shù)可視化，獲得空間中磁性體的概率分布。

如圖1所示，在笛卡爾坐標(biāo)系中，取平面z=0為觀測平面。假設(shè)在觀測平面上某一點（xi，yi，zi）測得的磁梯度張量為Bzz（xi，yi，zi）。觀測平面下的一點（ξ，η，ζ）存在一個小的長方體形的磁性體q，把它近似成一個磁偶極子，計算出它在點（xi，yi，zi）產(chǎn)生的磁異常為Bzz，q（xi，yi，zi）。

圖1 相關(guān)成像原理示意圖Fig.1 Schematic diagram of correlation imaging method

定義觀測平面內(nèi)測得的磁異常Bzz（xi，yi，zi）與磁性長方體q在觀測平面上產(chǎn)生的磁異常Bzz，q（xi，yi，zi）的皮氏（Pearson）互相關(guān)系數(shù)為

式中，N為測量點總數(shù)。

根據(jù)柯西不等式，可知Cq的取值范圍是[–1，1]。Cq表征了實測磁異常與磁性體q產(chǎn)生的磁異常的相關(guān)性，即實測磁異常有多大的可能性是由磁性體q產(chǎn)生的。Cq越接近于1，則該位置越有可能存在一個磁化率大的磁性體；Cq越接近于–1，則該位置越有可能不存在磁性體。

2 磁梯度數(shù)據(jù)的自適應(yīng)濾波

在二維測量平面上得到的磁梯度數(shù)據(jù)可以等效為一張二維圖像，通過快速傅里葉變換就能得到磁梯度數(shù)據(jù)的頻譜。

尺寸為M×N的函數(shù)f（x，y）的DFT為

式中，u和v為頻率變量。u=0，1，2，…，M–1；v=0，1，2，…，N–1。

給出F（u，v），可通過反DFT得到

式中，x和y為測量點的橫坐標(biāo)與縱坐標(biāo)。

圖2展示了不同深度的磁性體的磁梯度數(shù)據(jù)及幅度譜。在測量平面下方，淺層的小型磁性體引起的異常范圍小、變化快，對應(yīng)頻率域的高頻成分。深層較大的磁性體引起的異常范圍大、變化平緩，對應(yīng)頻率域的低頻成分。

圖2 不同深度的磁性體的磁梯度數(shù)據(jù)及幅度譜Fig.2 Magnetic gradient data and amplitude spectra of magnetic bodies at different depths

在磁梯度數(shù)據(jù)的相關(guān)成像過程中，如果成像空間中存在不同深度的多個磁性體，磁性體的磁場相互疊加，會導(dǎo)致成像結(jié)果的縱向分辨率嚴(yán)重降低，無法準(zhǔn)確估計磁性體的深度。但是由于不同深度磁性體的磁梯度數(shù)據(jù)在頻率域具有明顯的差異，因此可以利用頻率域濾波來分離不同深度的磁異常。

本文提出一種磁梯度數(shù)據(jù)的自適應(yīng)濾波方法，以改善相關(guān)成像的分辨率。具體步驟是：

1）測量磁性體的磁梯度張量數(shù)據(jù)Bzz（xi，yi）；

2）將測量平面以下的空間均勻剖分成M個網(wǎng)格單元；

3）用解析法計算第n個網(wǎng)格單元的磁梯度數(shù)據(jù)Bzz，q（xi，yi），并用快速傅里葉變換得到其頻譜Szz，q（ui，vi）；

4）將頻譜Szz，q（ui，vi）的每一個頻率分量按照幅度進行排序，保留幅度大小前50%的分量，標(biāo)記其余分量并置0，然后用逆傅里葉變換得到濾波后的磁梯度數(shù)據(jù)B’zz，q（xi，yi）；

5）濾除觀測數(shù)據(jù)Bzz（xi，yi）中被步驟4標(biāo)記的頻率分量，得到濾波后的觀測數(shù)據(jù)B’zz（xi，yi）；

6）計算B’zz，q（xi，yi）與B’zz（xi，yi）的相關(guān)系數(shù)；

7）重復(fù)步驟3–6依次計算M個網(wǎng)格單元的相關(guān)系數(shù)，得到三維相關(guān)成像結(jié)果。

該方法的特點是依據(jù)每一個網(wǎng)格單元的磁梯度數(shù)據(jù)頻譜來確定要保留測量數(shù)據(jù)中的哪些頻率成分，而不是采用固定的頻率閾值進行濾波。深處的磁性體產(chǎn)生的磁梯度張量低頻成分居多，就濾除觀測數(shù)據(jù)中的高頻成分，反之則濾除觀測數(shù)據(jù)中的低頻成分，然后再計算兩者的相關(guān)系數(shù)。

3 仿真模型實驗

本文設(shè)定反演空間的水平邊界范圍為–2.5～2.5 m，垂直邊界范圍為0～2.4 m。然后把反演空間劃分成立方體網(wǎng)格單元，每個立方體的大小都為0.05 m×0.05 m×0.05 m，共計480 000個網(wǎng)格單元。為確保數(shù)據(jù)能被完整觀測，觀測平面的范圍應(yīng)大于反演空間的水平邊界范圍。設(shè)定觀測平面的范圍為–4～4 m，均勻分布21×21個觀測點，觀測點間距為0.4 m。

為了驗證所提出的三維相關(guān)成像方法的有效性，在反演空間中設(shè)置了4個長方體模型來模擬水下磁性目標(biāo)。各個模型的參數(shù)設(shè)置如表1所示，模型的三維視圖如圖3所示。

表1 模型參數(shù)設(shè)置Table 1 Model parameter settings

圖3 模型的三維視圖Fig.3 3D view of model

以解析方法在觀測平面上進行正演，在高度為Z=0所計算的磁梯度張量分量理論值如圖4（a）所示?？紤]到實際測量環(huán)境中有噪聲的影響，在理論數(shù)據(jù)中加入2%的高斯白噪聲模擬含有磁測噪聲的信號，含噪聲信號如圖4（b）。

圖4 模型的磁梯度數(shù)據(jù)Fig.4 Magnetic gradient data of model

用快速傅里葉變換計算磁梯度數(shù)據(jù)的頻譜，結(jié)果如圖5所示。

圖5 磁梯度數(shù)據(jù)的幅度譜Fig.5 Amplitude spectrum of magnetic gradient data

采用第2章提出的自適應(yīng)濾波方法對圖4（b）所示的含噪聲信號進行相關(guān)成像，圖6顯示了成像結(jié)果。其中：圖6（a）與6（b）為Y=1.6 m處的剖面圖，圖6（c）與6（d）為X=1.4 m處的剖面圖，圖6（a）與6（c）為采用了所提濾波方法的結(jié)果，圖6（b）與6（d）為未采用濾波處理的結(jié)果。圖中黑色方框代表實際模型的位置，白色十字符號代表相關(guān)系數(shù)的極大值位置?？梢园l(fā)現(xiàn)，經(jīng)過濾波處理后，相關(guān)系數(shù)的極大值位置與模型的中心位置基本重合，只有圖6（c）中模型3在Z方向有0.05 m誤差；無濾波處理的成像結(jié)果中，如圖6（b），模型2在X和Z方向均有0.05 m誤差；圖6（d）中，模型3在Z方向均誤差0.1 m。

圖6 相關(guān)成像結(jié)果Fig.6 The correlation imaging results

由于實際應(yīng)用中需要通過相關(guān)系數(shù)的極大值位置估計磁性體中心的位置。用上述2種方法估計出4個模型的中心位置如表2所示。

表2 根據(jù)成像結(jié)果估計的模型中心位置Table 2 Estimated center position of model based on imaging results

由于相關(guān)成像方法需要將未知空間離散化為有限個網(wǎng)格，因此當(dāng)估計的模型中心位置與實際位置一致時，并不代表誤差為0。如表2中的模型1，估計的中心位置是（–1.6，–1.6，1.1），表明模型中心位于以（–1.6，–1.6，1.1）為中心，邊長為0.05 m的一個正方體網(wǎng)格內(nèi)，誤差區(qū)間為0～0.025 m。表2中各個模型的中心位置誤差范圍表明，提出的濾波方法對誤差有明顯改善，最多減小了1–0.25/1.06≈76%的誤差。

為了驗證低信噪比下所提濾波方法的效果，在理論信號中加入10%的高斯噪聲（如圖7所示）。由于每次實驗的高斯噪聲是隨機產(chǎn)生的，且10%的高斯噪聲隨機性已經(jīng)很大，分別在有無濾波處理的條件下進行了5次實驗，計算出磁性體位置誤差的平均值如表3所示。

表3 低信噪比下模型位置誤差Table 3 Model position errors under low SNR

圖7 含10%噪聲的磁梯度數(shù)據(jù)Fig.7 Magnetic gradient data with 10% noise

根據(jù)表3，濾波處理后模型1和4的位置誤差明顯小于無濾波處理的結(jié)果，而模型2和3的誤差沒有太大變化。原因可能是相關(guān)成像算法將磁性體等效為了球形的磁偶極子，模型2和3為長方體，用磁偶極子近似會有較大誤差，而模型1和4為正方體，更接近球形的磁偶極子。說明在低信噪比的情形下，濾波處理后對近似球形的模型定位結(jié)果更準(zhǔn)確。

在表1的基礎(chǔ)上將4個模型的深度增加0.6 m，并與圖4（b）所示信號一樣添加2%的高斯噪聲，對含噪聲信號進行相關(guān)成像以估計模型的中心位置，結(jié)果如表4所示。由于在成像結(jié)果X、Y、Z方向上均出現(xiàn)較多誤差，因此計算出3個方向上的誤差對結(jié)果進行評估。需要注意的是，表4中的模型在某一方向上的誤差值為0也不代表沒有誤差，而是理解成與表2一樣的誤差區(qū)間。

表4 模型加深后根據(jù)成像結(jié)果估計的模型中心位置Table 4 Center position of model estimated based on imaging results after deepening model

當(dāng)磁性體所處位置加深后，2種方法得到的結(jié)果均向上偏移，在Z方向上誤差較大，最高分別達到0.3 m和0.55 m，在水平方向上也都出現(xiàn)了0.05～0.1 m的誤差。說明在磁性體深度較深的情況下，本文所提濾波方法對成像效果的提升已不明顯，但相對未采用濾波處理的成像結(jié)果來說還是更接近真實情況。

4 結(jié)束語

針對磁梯度數(shù)據(jù)的相關(guān)成像方法縱向分辨率不足的問題，提出了一種自適應(yīng)濾波方法，根據(jù)不同深度磁性體的磁梯度數(shù)據(jù)頻譜差異使用不同的濾波閾值。仿真模型實驗表明，在平均深度1.1 m時，濾波處理后的成像結(jié)果能較為準(zhǔn)確反映磁性體的中心位置，有效減小了磁性體的定位誤差，最多可以減小76%；在平均深度1.7 m時，深度方向上最大誤差0.3 m，約為無濾波處理結(jié)果的54%；而且在低信噪比下，濾波方法還能減小近似球形的模型的定位誤差。總體而言，自適應(yīng)濾波可以提升相關(guān)成像法準(zhǔn)確定位磁性體的能力。