基于EMD的鐵磁性物質(zhì)檢測研究

董建榮，李丕丁，王 成

（上海理工大學(xué)健康科學(xué)與工程學(xué)院 上海 200093）

0 引言

核磁共振成像（Magnetic Resonance Imaging，MRI）作為臨床診斷的主要成像技術(shù)［1］，在工作過程中會產(chǎn)生強(qiáng)磁環(huán)境（0.2-3T），容易受鐵磁性物質(zhì)影響，甚至對被測試者身體安全產(chǎn)生危害。因此，必須對所有進(jìn)入MRI 工作環(huán)境中的人員進(jìn)行鐵磁性物質(zhì)檢測。

磁異常探測（Magnetic Anomaly Detection，MAD）作為一種被動探測技術(shù)，通過分析目標(biāo)物體磁性特征在環(huán)境磁場中產(chǎn)生的磁異信號檢測目標(biāo)物體。現(xiàn)已被廣泛應(yīng)用于航空探淺、水下鐵磁性物質(zhì)探測、地震預(yù)測、醫(yī)學(xué)核磁共振等領(lǐng)域［2］。目前，MAD 技術(shù)已與多種算法相結(jié)合，例如Pradhan 等［3］提出基于正交基函數(shù)的匹配濾波方法（Orthogonal Basis Functions，OBF），將磁異信號表示為3 個正交函數(shù)的線性組合匹配濾波。王磊等［4］提出基于自回歸（Auto Regression，AR）模型的白化濾波器，現(xiàn)將測量信號濾波處理，然后進(jìn)行OBF 檢波，最后利用最小熵檢測法（Minimum Entropy Detector，MED）和高階過零檢測法（High Order Croossing，HOC）檢測較低信噪比條件下的磁目標(biāo)。

然而，以上方法均建立在背景噪聲不變的前提下進(jìn)行實(shí)驗(yàn)，OBF 還需要獲取測量信號的先驗(yàn)信息。為解決以上問題，Nezhadshahbodaghi 等［5］提出基于能量檢測器的非抽樣離散小波變換法（Undecimated Discrete Wavelet Transform，UDWT）。該算法無需測量信號任何先驗(yàn)信息即可直接檢測磁異信號，但由于小波分析首先需要確定小波基，因此對不同特性目標(biāo)信號，可能顯示不同性能。此外，針對非線性、非平穩(wěn)性復(fù)雜信號，Inturi 等［6］提出經(jīng)驗(yàn)?zāi)B(tài)分解方法（Empirical Mode Decomposition，EMD），該方法將信號分解成多個固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF）以表征信號的基本信息，根據(jù)磁異信號特征自適應(yīng)調(diào)整分析精度，無需信號的任何先驗(yàn)信息和對背景噪聲進(jìn)行假設(shè)。

雖然MAD 方法已在各領(lǐng)域廣泛應(yīng)用，但在核磁共振成像領(lǐng)域鮮有研究。為此，本文采用MAD 方法檢測MRI工作環(huán)境中的鐵磁性物質(zhì)。在基于磁通門傳感器的鐵磁性探測系統(tǒng)基礎(chǔ)上，特征分析探測的鐵磁性物質(zhì)產(chǎn)生的磁異信號，并使用EMD算法從強(qiáng)噪聲背景中提取目標(biāo)磁異信號。

1 磁異常探測模型

本文設(shè)計磁探測系統(tǒng)處于靜止?fàn)顟B(tài)，檢測移動中的鐵磁性目標(biāo)物質(zhì)，假設(shè)鐵磁性目標(biāo)物體以恒定速度V沿直線運(yùn)動經(jīng)過磁探測系統(tǒng)，磁矩m在方向上和大小上恒定，CPA（closest proximity approach）為目標(biāo)物體距離磁探測系統(tǒng)的最小距離，當(dāng)目標(biāo)物體移動到CPA 位置附近，獲得的目標(biāo)信息最豐富。磁異常探測模型［7］如圖1所示。

Fig.1 Model of magnetic anomaly detection圖1 磁異常探測模型

由于地磁場作用，鐵磁性物質(zhì)會在區(qū)域范圍內(nèi)地磁場Be中產(chǎn)生擾動磁場Bi，從而產(chǎn)生磁場異常信號Ba，三者間的矢量關(guān)系如式（1）所示：

當(dāng)鐵磁性目標(biāo)物體直徑遠(yuǎn)小于其與磁探測系統(tǒng)的距離時，目標(biāo)磁場可等效于磁偶極子磁場［8］。根據(jù)畢奧薩伐爾定律，在遠(yuǎn)離磁矩m（目標(biāo)）處（距離為R）所產(chǎn)生的磁感應(yīng)強(qiáng)度Bm可表示為：

式中，m為偶極子磁矩，單位為A·m2，真空磁導(dǎo)率μ=4π× 10-7H/m，R為目標(biāo)物體與磁探測系統(tǒng)間的距離。因此，探測的磁異信號S可近似表示為磁偶極子場Bm在地磁場T（T>>Bm）上的投影［8］：

根據(jù)式（3），采用MATLAB 仿真軟件對磁異信號的等價模型進(jìn)行模擬?？紤]到地磁矢量在實(shí)際應(yīng)用中變化范圍有限，因此設(shè)為不變量［9］，磁異信號仿真波形如圖2所示。

Fig.2 Magnetic anomaly signal simulation waveform圖2 磁異信號仿真波形

當(dāng)檢測到磁異信號后，還需要能量檢測器探測磁異信號，對探測到的信號進(jìn)行平方求和，得到信號能量E，并將其與設(shè)定的能量閾值Em進(jìn)行比較，以完成磁異信號檢測［10］：

2.3 吸煙組、戒煙組及非吸煙組的椎體骨折陽性率比較 結(jié)果表明：吸煙組、戒煙組及非吸煙組的椎體骨折陽性率分別為23.3%(496/2 129)、21.8%(76/349)和14.2%(131/925)，組間差異有統(tǒng)計學(xué)意義(P<0.05)；經(jīng)兩兩比較，發(fā)現(xiàn)吸煙組與非吸煙組(P<0.001)、戒煙組與非吸煙組(P=0.001)的差異均有統(tǒng)計學(xué)意義。

由式（4）可知，當(dāng)E>Em時，表示探測系統(tǒng)檢測到磁異信號；當(dāng)E<=Em時，表示信號未達(dá)到設(shè)定能量閾值，無鐵磁性物質(zhì)磁場干擾。

2 EMD算法

經(jīng)驗(yàn)?zāi)B(tài)分解法能將任何復(fù)雜的時間、空間信號分解為從高頻到低頻的若干階固有模態(tài)函數(shù)（Intrinsic Mode function，IMF）［11］，每個IMF 均為穩(wěn)定的窄帶信號，分別表征信號的具體特征，而分解成固有模態(tài)函數(shù)的限制條件包括整個數(shù)據(jù)段極值個數(shù)和零點(diǎn)個數(shù)相差不超過1，在任何數(shù)據(jù)點(diǎn)上包絡(luò)（由局部極大值點(diǎn)形成）和下包絡(luò)（由局部極小值點(diǎn)形成）的平均值為0，即上下包絡(luò)局部對稱［12］。EMD 算法基本計算流程如圖3所示。

Fig.3 EMD algorithm flow chart圖3 EMD算法流程圖

由圖3 可見，EMD 算法的基本計算流程可分為以下4個步驟：

步驟1：取原始信號S(t)的局部極大與極小值點(diǎn)，分別采取三次樣條插值的方法構(gòu)建S(t)的上包絡(luò)Su(t)和下包絡(luò)Sl(t)。

步驟2：平均處理上包絡(luò)Su(t)和下包絡(luò)Sl(t)，獲得平均值函數(shù)m(t)，該值代表原始信號S(t)的低頻成分，通過S(t)減去m(t)得到新數(shù)據(jù)序列h1(t)。

步驟3：利用IMF 評判標(biāo)準(zhǔn)判斷h1(t)是否符合要求，若不符合，將h1(t)代入重復(fù)步驟1、步驟2，直至h1(t)達(dá)到IMF 標(biāo)準(zhǔn)，如此可得到表征S(t)最高頻率組成部分的第1個IMF 分量x1(t)。

步驟4：首先利用原始信號S(t)減去x1(t)得到剩余分量r1(t)。

最后，當(dāng)數(shù)據(jù)序列呈現(xiàn)單調(diào)特性或整個序列只存在一個極點(diǎn)的條件下，原始信號S(t)將終止分解。最終剩余分量即為EMD 分解所產(chǎn)生的殘差分量Res(t)，該分量代表原始信號的整體趨勢。綜上，EMD 將原始信號分解為若干個IMF 分量和殘差分量，即可用式（7）表示原始信號S(t)。

EMD 是一種時間尺度分析方法，對原始信號進(jìn)行EMD 分解所產(chǎn)生的若干IMF 分量，每個分量分別體現(xiàn)不同時間尺度特征，而該特征顯示的是非平穩(wěn)非線性信號的固有模態(tài)特征，因此需要選擇合適的分解層次進(jìn)行組合［13］，保證信號完整性，過濾高頻噪聲信號，獲取高信噪比目標(biāo)信號。本文在圖2 的磁異仿真信號上增加不同高斯噪聲，組成信噪比分別為-5dB、-7dB、-10dB 復(fù)雜背景下的磁異仿真信號，并使用EMD 算法對上述仿真信號進(jìn)行分解，獲取不同的IMF 分量并對其進(jìn)行篩選和重組［14］，由于高頻噪聲信號往往在高階IMF 分量中占比較高，因此在重組時將首先剔除高階IMF 分量。其中，一階、二階IMF 分量中高頻噪聲占比較大，故可選擇性地去除一、二、三階IMF 分量［15］，仿真結(jié)果如圖4所示。

圖4（a）、圖4（c）、圖4（e）分別為信噪比為-5dB、-7dB、-10dB 噪聲干擾的磁異仿真信號，經(jīng)過EMD 分解重構(gòu)的數(shù)據(jù)波形分別如圖4（b）、圖4（d）、圖4（f）所示。通過上述EMD 分解結(jié)果可知，對不同信噪比的磁異信號數(shù)據(jù)，EMD 均能較好地剔除高頻噪聲，提取磁異信號的準(zhǔn)確性較高，得到的處理信號具有較高的信噪比和穩(wěn)定性。但在實(shí)際探測環(huán)境中，由于MRI 受到強(qiáng)磁干擾、鐵磁性目標(biāo)物質(zhì)體積材質(zhì)不宜探測、人員體內(nèi)植入鐵磁性目標(biāo)物質(zhì)等因素，導(dǎo)致磁異信號被噪聲信號覆蓋。為了驗(yàn)證該算法的通用性和魯棒性，本文將采用實(shí)測磁異信號數(shù)據(jù)對EMD 算法進(jìn)行實(shí)驗(yàn)。

Fig.4 EMD decomposition and reconstruction results of different SNR simulation signals圖4 不同SNR仿真信號的EMD分解重構(gòu)結(jié)果

3 仿真實(shí)驗(yàn)

3.1 數(shù)據(jù)采集

本文實(shí)驗(yàn)數(shù)據(jù)由實(shí)驗(yàn)室自主設(shè)計的基于磁通門傳感器的鐵磁性物質(zhì)探測系統(tǒng)收集，該系統(tǒng)通過ADC 模塊將模擬磁異信號轉(zhuǎn)化為數(shù)字信號，配置了藍(lán)牙串口傳輸功能的上位機(jī)以顯示和保存數(shù)據(jù)。數(shù)據(jù)采集地點(diǎn)為不同醫(yī)院核磁共振室門外，以鉗子、手機(jī)、鑰匙和硬幣共4 種常見的鐵磁性物質(zhì)被測目標(biāo)。首先固定被測目標(biāo)與探測系統(tǒng)的距離，測試探測系統(tǒng)在不同角度下的探測性能；然后固定被測目標(biāo)與探測系統(tǒng)的角度，測試探測系統(tǒng)在不同距離的探測性能。

3.2 EMD算法驗(yàn)證

在測試實(shí)驗(yàn)中，為了避免復(fù)雜背景環(huán)境會給磁異信號與無關(guān)因素對檢測造成干擾，采取控制變量法進(jìn)行測試。同時，為避免實(shí)驗(yàn)數(shù)據(jù)量較大，本文僅以4 種被測物體在0°測試角度、120cm（鉗子、手機(jī)的最遠(yuǎn)測試距離）和50cm（鑰匙、硬幣的最遠(yuǎn)測試距離）距離條件下進(jìn)行測試，在MATLAB 中編寫EMD 數(shù)據(jù)處理的算法模型，使用UDWT算法對數(shù)據(jù)進(jìn)行濾波處理，處理結(jié)果如圖5所示。

Fig.5 EMD processing results of the magnetic abnormal signals of four measured substances圖5 4種被測物質(zhì)磁異信號的EMD處理結(jié)果

圖5（a）、圖5（d）、圖5（g）、圖5（j）分別為鉗子、手機(jī)、鑰匙、硬幣在MRI 室外復(fù)雜環(huán)境下由鐵磁性探測系統(tǒng)探測的磁異信號波形，可明顯看到信號存在較大噪聲干擾。由于被測物體對磁場的影響程度不同，因此造成平均幅值依次降低。圖5（b）、圖5（e）、圖5（h）、圖5（k）分別為上述磁異信號經(jīng)過EMD 處理過后重構(gòu)的信號波形，圖中圓圈部分波形為被測目標(biāo)所產(chǎn)生的磁異信號。圖5（c）、圖5（f）、圖5（i）、圖5（l）分別為上述磁異信號經(jīng)過UDWT 處理后的信號波形，經(jīng)過UDWT 濾波后雖保留了信號的有效信息，但信號仍然存在毛刺。

通過比較可知，EMD 處理后的波形相較于UDWT 濾波后更平滑，在去除高頻噪聲信號的同時保留了原始信號的有效成分，保證了原始信號的完整性。由測試結(jié)果可知，EMD 算法相較于傳統(tǒng)方法，處理不同測試物體的磁異信號更穩(wěn)定，能有效過濾提取MRI 復(fù)雜背景環(huán)境下鐵磁性物質(zhì)的磁異信號。此外，本文采用信號信噪比SNR 作為數(shù)據(jù)處理效果衡量指標(biāo)，并將EMD 方法與UDWT 方法所處理的結(jié)果進(jìn)行比較，結(jié)果如圖6-圖7所示。

Fig.6 SNR of EMD processing signals圖6 EMD處理信號的信噪比

Fig.7 SNR of UDWT processing signals圖7 UDWT處理信號的信噪比

圖6 中4 條折線分別代表4 種被測物體，橫坐標(biāo)為測試距離，縱坐標(biāo)為信噪比。由此可見，磁異信號經(jīng)過EMD處理過后的信噪比高達(dá)到16dB 以上，相較于原始信號具有較大提升。由圖7 可見，UDWT 方法對磁異信號的濾波效果不佳，濾波后信號的信噪比處于8dB 以下，且當(dāng)被測物體距離增加，磁異信號強(qiáng)度會隨之減小，信噪比處于下降趨勢。

綜上所述，EMD 方法相較于UDWT 方法，去除高頻噪聲更高效，能大幅度提升信號信噪比，對磁異信號的探測效果良好。同時，該方法可提高探測距離與靈敏度，更好的保證信號完整性，為檢測鐵磁性目標(biāo)物質(zhì)提供良好的數(shù)據(jù)保證。

4 EMD方法改進(jìn)

目前，EMD 方法存在模態(tài)混疊問題。模態(tài)混疊是指同一IMF 分量包含差異較大的特征時間尺度，或相近特征時間尺度分布在不同的IMF 分量中，該問題主要由于在分解的過程中局部極值在短時間內(nèi)發(fā)生多次跳變所導(dǎo)致，將致使IMF 分量失去物理意義，影響后續(xù)時頻分析處理。

為此，鄒瑛珂等［16］提出協(xié)方差疊加經(jīng)驗(yàn)?zāi)B(tài)分解（Covariance Superposition Empirical Mode Decomposition，CSEMD）的改進(jìn)型EMD 方法，基本原理是加強(qiáng)信號中具有高相關(guān)性的信號，減少低相關(guān)性的信號（一般為噪音）。CSEMD 方法首先對原信號進(jìn)行EMD 分解取得IMF 分量；然后對各IMF 分量和原信號進(jìn)行協(xié)方差運(yùn)算，將協(xié)方差標(biāo)準(zhǔn)化后進(jìn)行排序，將最大值與最小值相減后取平均；最后在原信號基礎(chǔ)上增加平均值進(jìn)行EMD 分解得到高信噪比的IMF 分量。

本文實(shí)驗(yàn)測試數(shù)據(jù)共200 組，4 種被測鐵磁性物質(zhì)各50 組，由于磁異信號與探測角度、探測距離、被測物質(zhì)的鐵磁性成分占比相關(guān)，因此固定測量角度，依次增加測試距離進(jìn)行實(shí)驗(yàn)。其中，鉗子與手機(jī)的測試距離為10cm、30cm、50cm、80cm、120cm，鑰匙與硬幣的測試距離為10cm、20cm、30cm、40cm、50cm，每個均距離測試10次。

首先將實(shí)驗(yàn)數(shù)據(jù)分別進(jìn)行CSEMD、EMD、UDWT 處理，然后通過能量檢測器檢測磁異信號。若無磁異信號時，能量檢測器的輸出較為平緩；當(dāng)出現(xiàn)磁異信號時，檢測器表現(xiàn)為能量波形凸起；若能量幅值超過設(shè)定能量閾值，則將觸發(fā)光聲報警。鐵磁性物質(zhì)檢出率如表1所示。

Table 1 Detection rate of different ferromagnetic substances by three methods表1 3種方法對不同鐵磁性物質(zhì)的檢測率 （%）

由表1 可知，UDWT 方法對于磁異信號的檢測效率較低，并且隨著鐵磁性物質(zhì)含鐵磁量減小而降低，無法有效檢測被測鐵磁性物體；EMD 方法在磁異信號的檢測方面明顯優(yōu)于UDWT 方法，大幅度提升了檢測鐵磁性物質(zhì)的準(zhǔn)確率；CSEMD 方法在EMD 的基礎(chǔ)上對模態(tài)混疊問題有所抑制，提高了磁異信號的檢測效率。

5 結(jié)語

本文提出無需任何先驗(yàn)信息且適用于非平穩(wěn)、非線性信號的EMD 算法，并將該算法應(yīng)用于MRI 室外強(qiáng)磁環(huán)境中檢測微弱的磁異信號。首先對包含高頻噪聲的磁異信號進(jìn)行EMD 分解，剔除含有高頻噪聲的IMF 分量再重構(gòu)其他分量以保證信號的完整性；然后使用改進(jìn)型CSEMD 方法抑制EMD 方法模態(tài)混疊問題；最后通過能量檢測器檢測磁異信號，提高EMD 方法對鐵磁性物體檢測準(zhǔn)確率。

實(shí)驗(yàn)數(shù)據(jù)表明，EMD 算法將磁異信號波形的信噪比提高至16dB 以上，相較于傳統(tǒng)的UDWT 算法提高了8dB 以上。同時，使用改進(jìn)EMD 方法檢測磁性物質(zhì)的準(zhǔn)確率高達(dá)90%以上，能有效檢測MRI 復(fù)雜環(huán)境背景下微弱的磁異信號，顯著提高了鐵磁性探測系統(tǒng)的檢測距離和檢測效率。然而，本文僅將該方法運(yùn)用于MRI 領(lǐng)域，下一步將在不同領(lǐng)域測試算法的適用性。