董建榮，李丕丁，王 成

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

0 引言

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

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

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

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

1 磁異常探測(cè)模型

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

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

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

當(dāng)鐵磁性目標(biāo)物體直徑遠(yuǎn)小于其與磁探測(cè)系統(tǒng)的距離時(shí)，目標(biāo)磁場(chǎng)可等效于磁偶極子磁場(chǎng)［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)物體與磁探測(cè)系統(tǒng)間的距離。因此，探測(cè)的磁異信號(hào)S可近似表示為磁偶極子場(chǎng)Bm在地磁場(chǎng)T（T>>Bm）上的投影［8］：

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

Fig.2 Magnetic anomaly signal simulation waveform圖2 磁異信號(hào)仿真波形

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

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

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

2 EMD算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 EMD算法驗(yàn)證

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

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

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

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

Fig.6 SNR of EMD processing signals圖6 EMD處理信號(hào)的信噪比

Fig.7 SNR of UDWT processing signals圖7 UDWT處理信號(hào)的信噪比

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

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

4 EMD方法改進(jìn)

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

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

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

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

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

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

5 結(jié)語

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

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