基于總體最小二乘的核四極矩共振參數(shù)估計(jì)

朱凱然 何學(xué)輝 鄭小保 蘇 濤

（1.西安電子科技大學(xué)雷達(dá)信號處理國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710071；2.中國電波傳播研究所，山東 青島 266107）

引 言

進(jìn)入21世紀(jì)后，反恐安全檢測技術(shù)已成為人們的熱點(diǎn)話題，同時(shí)也成為國防安全中重要的一部分，其目的是要把恐怖活動遏制在未遂狀態(tài)?？v觀這些恐怖活動，爆炸仍舊是恐怖分子使用的最主要方式，對爆炸物的及時(shí)有效探測是反恐活動的重要環(huán)節(jié)。

核四極矩共振 （NQR）是一種固態(tài)射頻技術(shù)［1-9，11-16］，它可以用于檢測四極矩原子核的存在性，例如普遍存在于爆炸物中的氮14原子核，它可作為一種判決依據(jù)對爆炸物進(jìn)行探測。通過檢測氮14原子核的NQR信號，不但可以判定爆炸物的存在性，而且還可以根據(jù)NQR譜線特征確定是何種爆炸物，可以說NQR信號是爆炸物的“指紋”。

NQR信號通常分為兩種類型。一種稱為自由感應(yīng)衰減（FID）信號［12］，它是由單脈沖激勵得到的信號；另一種稱為自旋回波串信號［6］，它是由激勵脈沖串激勵得到的信號，可以在短時(shí)間內(nèi)使NQR信號獲得更高的信噪比。然而，NQR信號本身非常弱，并且易受線圈內(nèi)熱噪聲和外部射頻干擾（RFI）的影響，故此NQR信號極低信噪比制約了NQR技術(shù)在爆炸物探測方面的應(yīng)用。文獻(xiàn)［1］應(yīng)用非線性最小二乘法提高信噪比，該方法在已知信源個(gè)數(shù)的前提下，建立NQR信號模型，在頻率維和衰減系數(shù)維進(jìn)行窮舉搜索，計(jì)算量大，且需要先驗(yàn)信息。文獻(xiàn)［9］運(yùn)用基于最小二乘法的線性預(yù)測方程估計(jì)NQR信號的參數(shù)。本文采用基于總體最小二乘法［10］的線性預(yù)測方程對NQR信號的參數(shù)進(jìn)行估計(jì)，該方法大大降低了計(jì)算量，且無需先驗(yàn)信息；在計(jì)算線性預(yù)測系數(shù)時(shí)，它克服了普通最小二乘法未考慮擾動對數(shù)據(jù)矩陣的影響，獲得了更好的估計(jì)性能。

1.信號模型

窄帶接收機(jī)接收到的FID信號可以表示成一個(gè)指數(shù)衰減的正弦信號［1］，表達(dá)式為

式中:A、β、ω和θ分別為FID信號的幅度、衰減系數(shù)、頻率以及初相；ω通常是環(huán)境參數(shù)的函數(shù)，如溫度。假設(shè)

式中:f，fs分別表示FID信號頻率和采樣率。式（1）變?yōu)?/p>

實(shí)際觀測數(shù)據(jù)可以表示為

式中:IN，wN分別表示干擾信號和噪聲信號，均為N×1的向量；（·）T表示轉(zhuǎn)置。式（4）可以看作在色噪聲背景下，F(xiàn)ID信號參數(shù)估計(jì)問題。

2.基于總體最小二乘的參數(shù)估計(jì)

在求解矩陣方程Ax=b的時(shí)候，普通最小二乘是用一個(gè)范數(shù)平方為最小的擾動向量e去干擾數(shù)據(jù)向量b，以校正b中存在的噪聲，當(dāng)A和b二者均存在干擾時(shí)，求解矩陣方程Ax=b的最小二乘方法將會導(dǎo)致大的方差，此時(shí)，需要同時(shí)考慮矩陣A和向量b中的擾動，即求解矩陣方程（A+E）x=b+e，這就是總體最小二乘方法的基本思想。由式（4）構(gòu)造線性預(yù)測方程為

A為（N-L）×L 的 Hankel矩陣；x=（x1，x2，…，xL）T為線性預(yù)測向量；b=（y0，y1，…，yN-L-1）T.構(gòu)造增廣矩陣B=［b，A］，對矩陣B進(jìn)行奇異值分解（SVD）

式中:U=［u1，u2，…，uN-L］∈C（N-L）×（N-L）；均為正交矩陣；奇異值矩陣

其奇異值按照降序排列，即σ1≥σ2≥…≥σq＞σq+1≈…≈σL+1，且q≤min（N-L，L+1）是增廣矩陣B的有效秩；（·）H為共軛轉(zhuǎn)置。那么，式（5）的總體最小二乘解為

式中:S-（q）是矩陣的S（q）逆矩陣，且

式中:S（q）為（L+1）×（L+1）維矩陣；νi是酉矩陣V的第i列。

得到xTLS后，由式（5）中任意一個(gè)差分方程可得特征多項(xiàng)式

解得式（11）的根zi，i=1，…，L，由FID信號的衰減特性知，尋找 Re（ln（zi））＜0，i=1，…，L 的根zi，i=1，…，p，那么，式（3）可簡化為幅度參數(shù)a的線性方程，用矩陣形式表示為

式中，Z是N×p維Vandermonde矩陣。由于zi各不相同，那么Vandermonde矩陣Z的各列線性獨(dú)立，可得式（12）的最小二乘解為

式中:Z?=（ZHZ）-1ZH是Z的偽逆矩陣。依據(jù)

可求得幅值、相位、衰減因子和頻率。

得到p組信號分量估計(jì)參數(shù)后，依據(jù)最大似然估計(jì)準(zhǔn)則

式中:Zi是Z的第i列；‖·‖2為2-范數(shù)。對應(yīng)于J最小的那組參數(shù)，便是所求參數(shù)。式（17）的最小值可以轉(zhuǎn)化為求式（18）的最大值，即

式中:Πz=Z（ZHZ）-1ZH.

此時(shí)，可以獲得未知噪聲模型下的信噪比估計(jì)為

SNR′或稱之為檢測增益，它表征了NQR信號能量與背景噪聲能量之比。

值得注意的是，式（5）的最佳線性預(yù)測器的參數(shù)等效于自回歸（AR）模型的參數(shù)，最小預(yù)測誤差功率也就等于激勵的白噪聲功率，此時(shí)的線性預(yù)測器也有白化濾波的作用。式（9）中減掉Δσ項(xiàng)，也具有消除噪聲的作用。本文方法的步驟可歸納如下:

步驟1 根據(jù)式（5）與式（6）構(gòu)造線性預(yù)測方程，線性預(yù)測器階數(shù)通常取L≈N／3.

步驟2 計(jì)算增廣矩陣B的SVD，并存儲右奇異矩陣V；確定B的有效秩q；利用式（10）和式（9）計(jì)算矩陣S（q）；求S（q）的逆矩陣S-（q）；由式（8）得到線性預(yù)測器的系數(shù)xTLS.

步驟3 由式（11）解得特征多項(xiàng)式的根，zi，i=1，…，L，依據(jù) Re（ln（zi））＜0，i=1，…，L 可以得到根zi，i=1，…，p，構(gòu)造N×p維Vandermonde矩陣Z.

步驟4 由式（15）得到幅度信息的最小二乘解。

步驟5 由式（16）可以得到幅值、相位、衰減因子以及頻率信息，依據(jù)式（18）、式（19）分別可以確定FID信號的參數(shù)和檢測增益。

3.計(jì)算機(jī)仿真結(jié)果

應(yīng)用基于總體最小二乘法的線性預(yù)測方法估計(jì)得到的NQR參數(shù)結(jié)果與文獻(xiàn)［9］的方法進(jìn)行比較。從仿真數(shù)據(jù)和實(shí)測數(shù)據(jù)兩方面進(jìn)行仿真，并比較上述兩種方法的參數(shù)估計(jì)性能。

3.1 仿真數(shù)據(jù)

在室溫下，采集背景噪聲數(shù)據(jù)，然后加入相對于背景噪聲不同信噪比（SNR）下的NQR模擬信號，生成仿真數(shù)據(jù)。根據(jù)黑索金（RDX）實(shí)測數(shù)據(jù)給出NQR模擬信號的參數(shù):頻率為2.3kHz，衰減常數(shù)為0.0064，初相位為30°，采樣率為156.25kHz，信噪比從-5dB變化到20dB，步長為0.5dB.每個(gè)信噪比對背景噪聲樣本（22000個(gè)脈沖）進(jìn)行100次隨機(jī)不重復(fù)抽取（2000個(gè)脈沖）累加實(shí)驗(yàn)。

隨著信噪比的變化，參數(shù)估計(jì)值如圖1所示。其中，圖1（a）為幅度信息估計(jì)，在高信噪比下，兩者估計(jì)值均接近于理論值；圖1（b）為衰減系數(shù)估計(jì)，隨著信噪比的增加，估計(jì)值趨于真實(shí)值，若以衰減系數(shù)為檢測變量，那么在較低信噪比下，本文方法優(yōu)于文獻(xiàn)［9］的方法；圖1（c）為頻率估計(jì)值，隨著信噪比的增加，估計(jì)值也趨于真實(shí)值，本文方法優(yōu)于文獻(xiàn)［9］的方法；圖1（d）為初始相位估計(jì)值，隨著信噪比的增加，估計(jì)值也趨于真實(shí)值；文獻(xiàn)［9］方法估計(jì)值更接近于真實(shí)值。圖1（e）為檢測增益，本文方法得到的檢測增益優(yōu)于文獻(xiàn)［9］方法約0.8dB.

3.2 實(shí)測數(shù)據(jù)

實(shí)測數(shù)據(jù)源于研究RDX的FID信號特性實(shí)驗(yàn)［13-16］。探測RDX譜線頻率為3.411MHz，在該共振頻率下，其自旋-晶格弛豫時(shí)間在環(huán)境溫度為283～295K 條件下大約為12ms［6，8］，自旋-自旋弛豫時(shí)間約為8ms.室溫下，將300g RDX置于封閉的感應(yīng)線圈內(nèi)，用激勵脈沖寬度為130μs，脈沖間隔為64.1ms（約大于5倍自旋-晶格時(shí)間）的單脈沖序列激勵RDX.中心頻率為100kHz，采樣率為1.25 MHz.

處理中，選取第一個(gè)數(shù)據(jù)文件內(nèi)第一個(gè)脈沖為起始脈沖，以步長為62個(gè)脈沖順序增加脈沖數(shù)，進(jìn)行脈沖相干積累。理論上，N個(gè)等幅相干信號進(jìn)行相干積累，可以使輸出信噪比改善N倍。但實(shí)際上，由于相鄰回波信號的相位相干性難以保證，因此信噪比改善一般達(dá)不到N倍［1，11-12］。不同積累脈沖數(shù)下，兩種方法的FID信號參數(shù)估計(jì)性能比較如圖2所示。圖2（a）和圖2（d）說明了脈沖相干積累的可能性，因?yàn)殡S著脈沖數(shù)的增加，幅度信息也隨之線性增加，初始相位趨于穩(wěn)定，初始相位的平均值為35.82°.圖2（b）為FID信號的衰減系數(shù)估計(jì)，趨于0.007，而文獻(xiàn)［9］方法表現(xiàn)出較大的波動。圖2（c）為FID信號的頻率估計(jì)，趨于2.35kHz.圖2（e）為信噪比估計(jì)值，當(dāng)積累脈沖數(shù)小于2500個(gè)時(shí)，對應(yīng)的信噪比估計(jì)值小于10dB，文獻(xiàn)［9］方法低于本文方法，大于2500個(gè)時(shí)，兩種方法的性能一致，且趨于線性增加。綜上，在較低信噪比下，與基于普通最小二乘法的文獻(xiàn)［9］方法相比，本文方法呈現(xiàn)出較好的參數(shù)估計(jì)性能，同時(shí)也驗(yàn)證了FID信號模型的正確性。

4.結(jié) 論

在研究RDX的FID信號特性的基礎(chǔ)上，針對NQR信號參數(shù)的估計(jì)問題，應(yīng)用基于總體最小二乘法的線性預(yù)測方法，并考慮了線性預(yù)測方程左右擾動矩陣的影響，獲得了較好的參數(shù)估計(jì)性能。但該類線性預(yù)測方法也有不足之處，當(dāng)數(shù)據(jù)長度較大時(shí)，數(shù)據(jù)矩陣的行數(shù)太大，且為了得到好的預(yù)測系數(shù)，通常需要很大的線性預(yù)測階數(shù)。NQR技術(shù)是一種有效的檢測爆炸物技術(shù)，具有良好的發(fā)展前景，NQR信號極低信噪比一直是制約該技術(shù)在爆炸物探測方面應(yīng)用的瓶頸，研究低信噪比下爆炸物檢測算法成為亟待解決的問題。下一步的工作，將進(jìn)一步研究FID信號特性，穩(wěn)態(tài)自由振蕩脈沖序列（SSFP）回波信號特性，自旋鎖定脈沖序列（SLMP）回波信號特性以及低信噪比下爆炸物檢測算法。

［1］JAKOBSSON A，MOSSBERG M，ROWE M D，et al.Exploiting temperature dependency in the detection of NQR signals［J］.IEEE Transactions on Signal Processing，2006，54（5）:1610-1616.

［2］TAN Yingyi，TANTUM S L and COLLIN L M.Landmine detection with nuclear quadrupole resonance ［C］／／IEEE International Geoscience and Remote Sensing Symposium，2002，3:1575-1578.

［3］SOMASUNDARAM S D，JAKOBSSON A，and BUTT N R.Countering radio frequency interference in single-sensor quadrupole resonance［J］.IEEE Geoscience And Remote Sensing Letters，2009，6（1）:62-66.

［4］SOMASUNDARAM S D，JAKOBSSON A，ROWE M D，et al.Robust detection of stochastic nuclear quadrupole resonance signals［J］.IEEE Transactions on Signal Processing，2008，56（9）:4221-4229.

［5］TAN Yingyi，TANTUM S L，and COLLIN L M.Cramer-Rao lower bound for estimating quadrupole resonance signals in Non-Gaussian noise［J］.IEEE Signal Processing Letters，2004，11（5）:490-493.

［6］RUDAKOV T N，MIKHALTSEVITCH V T，and FLEXMAN J H.Modified steady-state free precession pulse sequences for the detection of pure nuclear quadrupole resonance［J］.Solid State Nuclear Magnetic Resonance，2004，25（1-3）:94-98.

［7］GARROWAY A N，BUESS M L，MILLER J B，et al.Remote sensing by nuclear quadrupole resonance［J］.IEEE Transactions on Geoscience and Remote Sensing，2001，39（6）:1108-1118.

［8］SMITH J A S.Nitrogen-14quadrupole resonance detection of RDX and HMX based explosives［C］／／European Convention on Security and Detection.London，1995:288-292.

［9］GRECHISHKIN V S，GRECHISHKINA R V，and HEO H.Signal prcessing methods in NQR ［C］／／Explosives Detection Using Magnetic and Nuclear Resonance Techniques，NATO Science for Peace and Security Series B:Physics and Biophysics.2009:159-170.

［10］張賢達(dá).矩陣分析與應(yīng)用［M］.北京:清華大學(xué)出版社，2004:403-452.

［11］趙振維，婁 揚(yáng)，金燕波，等.基于自適應(yīng)濾波技術(shù)的 NQR信號處理［J］.電波科學(xué)學(xué)報(bào)，2008，23（3）:429-433.ZHAO Zhenwei，LOU Yang，JIN Yanbo，et al.Signal processing for NQR based on adaptive filtering［J］.Chinese Journal of Radio Science，2008，23（3）:429-433.（in Chinese）

［12］原 普，毛云志，郭華民.爆炸物檢測中微弱NQR信號的處理［J］.電波科學(xué)學(xué)報(bào)，2006，21（2）:204-208.YUAN Pu，MAO Yunzhi，GUO Huamin.Weak signal receiving and processing for nuclear quadrupole resonance detection［J］.Chinese Journal of Radio Science，2006，21（2）:204-208.（in Chinese）

［13］夏佑林，葉朝輝.14N核四極共振中自旋鎖定的偏共振效應(yīng)［J］.物理學(xué)報(bào)，1995，44（6）:970-976.XIA Youlin，YE Chaohui.Offset effects of spin-locking in14N nuclear quadrupolar resonance［J］.Acta Physica Sinica，1995，44（6）:970-976.（in Chinese）

［14］夏佑林，葉朝輝.14N核四極共振自旋系統(tǒng)自旋-晶格弛豫時(shí)間的測量［J］.波譜學(xué)雜志，1994，11（4）:327-355.XIA Youlin，YE Chaohui.The measurments of spinlattice relaxation times for a14N nuclear quadrupolar resonance spin system［J］.Chinese Journal of Magnetic Resonance，1994，11（4）:327-355.（in Chinese）

［15］夏佑林，葉朝輝.自旋I=1和3／2的粉末樣品的核四極共振波譜學(xué)—I.對脈沖的響應(yīng)［J］.自然科學(xué)進(jìn)展，1995，5（4）:467-473.

［16］夏佑林，葉朝輝.自旋I=1和3／2的粉末樣品的核四極共振波譜學(xué)—II.核四極共振回波及應(yīng)用［J］.自然科學(xué)進(jìn)展，1995，5（5）:559-464.