姬婷婷， 張 杰， 王國宇

(中國海洋大學信息科學與工程學院，山東 青島 266100)

矢量序列信號的瞬時頻率分析具有重要的應(yīng)用需求，諸如雷達海雜波分析、海洋聲信號分析等。瞬時頻率計算對噪聲污染非常敏感，因此矢量序列噪聲抑制是不可回避的第一步環(huán)節(jié)。但是，從噪聲污染信號中提取信號，面臨的是一個“病態(tài)”投影信號復(fù)原的挑戰(zhàn)性問題。對于海雜波尤其如此。噪聲污染又加上海面結(jié)構(gòu)所導(dǎo)致的信號非線性和非平穩(wěn)性，更加劇了信號分析的困難。

現(xiàn)代信號處理越來越依賴于概率和統(tǒng)計以解決信號處理中的挑戰(zhàn)問題[1]。對于矢量形式的高斯噪聲，人們認識到它的相位服從均勻分布(白色)而矢徑服從高斯分布。借助信息熵的術(shù)語，均勻分布的信息熵最大，也意味著它的隨機性最強，因此需要最為精細尺度的處理方法。而對于服從非均勻分布的矢徑序列，有效的途徑是處理算法的尺度能夠依賴于數(shù)據(jù)的局部性尺度。

縱觀近年來已有的分解處理方法，全變分TV(Total Variation)方法具有最為精細的時間變化分析性能。全變分降噪模型[2]自從提出之后，由于物理意義的明顯和精細處理特點，在圖像分解、復(fù)原、分割等多方面得到廣泛應(yīng)用，目前已經(jīng)形成了一種信號和圖像處理領(lǐng)域的藝術(shù)級別的流派。TV方法的數(shù)字化計算過程中有2個參數(shù)，調(diào)整因子和迭代次數(shù)，需要進行選擇，文獻[2]給出了這參數(shù)調(diào)整因子的計算公式。而后，諸多學者研究了對處理效果的影響、諸如圖像的邊緣保持特性等。文獻[3]在超聲信號降噪的TV應(yīng)用中認為該參數(shù)計算的耗費太大，經(jīng)驗性的給出了該參數(shù)的選擇范圍。

1998年，Nudun Huang等[4]提出了適應(yīng)于非平穩(wěn)非線性信號的EMD(Empirical Mode Decomposition)分解方法。EMD將信號看作為“零均值快變震蕩信號與慢變震蕩信號的疊加”，其分解的主要特點是依據(jù)待分解數(shù)據(jù)的“局部尺度”，將原始信號分解為不同震蕩頻率的細節(jié)和內(nèi)稟性模式IMFS。而后，通過Hilbert-Huang Spectrum計算信號的瞬時頻率。文中也證明了利用分解量可以完全重構(gòu)原始信號。但是，Nudun Huang的EMD方法僅僅適用于1D實數(shù)序列分解。

為了能夠進行矢量序列分析，近年來逐步發(fā)展了RIEMD[5]、BEMD[6]等基于EMD的分解方法。BEMD認為當分析數(shù)據(jù)是兩變量數(shù)據(jù)時，振動的表示是模糊的。他們以旋轉(zhuǎn)表示兩變量信號：“兩變量信號=快速旋轉(zhuǎn)信號疊加于慢速旋轉(zhuǎn)信號”之上。文獻[6]同時給出了包絡(luò)均值的求取算法。文獻[7]注意到，對于規(guī)定兩個方向的情況RIEMD與BEMD 是等價的。文獻[8]應(yīng)用BEMD研究了在認知雷達場景分析的應(yīng)用。

根據(jù)噪聲統(tǒng)計特性及TV和EMD分析方法的特點，本文提出由TV相位降噪和矢徑EMD分解相組合的矢量序列降噪、分解方法。在相位降噪中，以最大信雜比為準則，優(yōu)化選擇迭代計算中的調(diào)整因子和迭代次數(shù)，以獲得最大信雜比的降噪效果。在矢徑EMD分解中，為了防止重構(gòu)矢量的相位跳變，對分解增加了非負判定環(huán)節(jié)。最后，將經(jīng)過降噪的相位序列和EMD分解的矢徑序列，以采樣時間(樣點序號)為參考，一一對應(yīng)組合，重構(gòu)矢量。考慮到歷史的傳承性，我們將該種矢量分解方法稱之為VEMD(Vector Empirical Mode Decomposition)方法。

1 矢量時間序列的分解

1.1 相位序列TV降噪

任意一個矢量時間序列可以表示為

Z(t)=a(t)exp{jφ(t)}。

(1)

其中：a(t)是矢量的矢徑；φ(t)是矢量在時刻t的相位。

將噪聲污染的相位序列φ(t)表達為

φ=y+n。

(2)

式中n為噪聲。TV 降噪模型為

(3)

式中Ω信號定義域，調(diào)整參數(shù)λ>0控制平滑量。(3)式中的第1項稱為正則項，第2項稱為保真項。(3)式最速下降法的Euler-Lagrange方程是

(4)

(5)

TV降噪的計算過程涉及調(diào)整因子和迭代次數(shù)的選擇。以自適應(yīng)方法選擇將會導(dǎo)致大計算量耗費。最大信雜比是信號處理所追求的一般性準則。本文采用最大信雜比為準則的方法指導(dǎo)和迭代次數(shù)的選擇。信雜比估計采用了經(jīng)驗性估計方法。假設(shè)信號頻率低于噪聲頻率，將低頻功率譜峰面積與背景雜波的功率譜面積之比作為估計的信雜比。

最大信雜比TV相位降噪及參數(shù)過程如下：首先設(shè)定λ和的初始數(shù)值及其變化步長；在每一次迭代計算之后，估計信雜比。最后，以最大信雜比下的調(diào)整因子和迭代次數(shù)作為參數(shù)，應(yīng)用于TV降噪。一個實際過程的信雜比估計的圖形如1所示。

圖1 TV迭代過程參數(shù)變化的信雜比Fig.1 The SNR with the parameters adjusting in TV iteration processes

1.2 修改的EMD矢徑分解

考慮到矢量重構(gòu)，關(guān)于矢徑分解的一個重要約束是不能夠改變原數(shù)值的符號(正/負)。因為改變一個矢量的正負，等價于矢量相位的變化。本文在原始EMD分解框架的基礎(chǔ)上，增加一個IMF的判別環(huán)節(jié)，?lmf(ti)<0，則分解停止。

EMD分解框架：

(2)求局部包絡(luò)：在相鄰的極大值之間(極小值之間同樣處理)，利用立方樣條插值，形成極大值包絡(luò)emax(和極小值包絡(luò)emin)。

(3)局部均值曲線：對每一個樣點，均值曲線：

(6)

(4)Imf非負判斷：如果?Imf(ti)<0, 停止。

(5)從信號中減均值，得到分解的細節(jié)Detail(i)=x(i)-Imf(i)。

經(jīng)過上述過程，從原始信號得到第一層分解的Imf1和Detail1。

將Imf1作為新的待分解信號，重復(fù)(1)～(5)，直到L層的ImfL。其中，(4)條為本文加入的約束條件。

1.3 矢量重構(gòu)

將VEMD分解的最后一層內(nèi)稟性模式與降噪相位序列φDenoise一一對應(yīng)，構(gòu)成矢徑、相位對，則稱為重構(gòu)的分解矢量序列，

(7)

它保留了原始信號的最大能量并且滿足TV相位保真性約束條件，因此將其看作為純凈信號的近似。

1.4 瞬時相位梯度與瞬時功率譜計算

(8)

為了計算瞬時功率譜，在(8)式的基礎(chǔ)上，將

(9)

賦給聯(lián)合坐標，作為瞬時功率譜。

2 分解實驗

為了驗證本文所提出方法的有效性，我們以人工產(chǎn)生的線性頻率調(diào)制信號，線性+正弦頻率調(diào)制的信號進行了分解噪聲實驗。實驗過程如下：首先生成純凈的矢量時間序列信號Zpure(t)，然后利用MATLAB功能函數(shù)，在信號上疊加信噪比為1dB的高斯白噪聲(請注意這是信雜比較低的情況)，然后利用本文的VEMD分解方法進行分解。為了表示的統(tǒng)一性，記矢徑分解的最后一層為Imf(t)。將Imf(t)與經(jīng)過TV降噪的相位序列一一對應(yīng)組合，構(gòu)成近似矢量Zre(t)=Imf(t)exp{jφDenoise(t)}，計算瞬時頻率。

實驗用例1，線性頻率調(diào)制信號Zpure(t)=3.0exp{j2π300(t2)/2}，實驗結(jié)果如圖2所示。

實驗用例2，線性頻率調(diào)制+正弦頻率調(diào)制信號，過程同實驗用例1。

Z(t)=3exp{j2π300(t2)/2+j5sin(2π10t)}實驗結(jié)果如圖3所示。

((a)純凈信號的瞬時頻率；(b)噪聲污染信號的瞬時頻率；(c)本文方法降噪信號的瞬時頻率；(d)純凈信號的瞬時功率；(e)噪聲污染信號的瞬時功率；(f)本文方法降噪信號的瞬時功率。(a) Instantaneous frequency of the pure signal ；(b) Instantaneous frequency of the noising signal；(c) The de-noising signal’s instantaneous frequency by method this paper；(d) Instantaneous power of the pure signal；(e) Instantaneous power of the noising signal；(f) Instantaneous power of the de-noising signal by method this paper.)

圖2 分解降噪實驗

Fig.2 Experiment de-noising by decomposition

由上述2個實驗結(jié)果的可見，本文的方法可以有效抑制矢量高斯噪聲。雖然在瞬時頻率的初始和最后階段有所形變，但是基本上保持了純凈信號的變化規(guī)律。

實驗用例3，加拿大IPIX雷達數(shù)據(jù)這里采用的試驗數(shù)據(jù)來自于加拿大McMaster大學的IPIX雷達數(shù)據(jù)庫網(wǎng)站，其雷達數(shù)據(jù)包括了HH、HV、VH，VV四種極化方式，每一個數(shù)據(jù)文件含有14個距離門的回波信號，其中目標為一個直徑1 m的球形密封器件，表面包裹了鋁箔，以增強反射回波信號。

((a)純凈信號的瞬時頻率；(b)噪聲污染信號的瞬時頻率；(c)本文方法降噪信號的瞬時頻率；(d)純凈信號的瞬時功率；(e)噪聲污染信號的瞬時功率；(f)本文方法降噪信號的瞬時功率。(a) Instantaneous frequency of the pure signal;(b) Instantaneous frequency of the noising signal;(c) The de-noising signal’s instantaneous frequency by method this paper;(d) Instantaneous Power of the pure signal;(e) Instantaneous Power of the noising signal;(f) Instantaneous Power of the de-noising signal by method this paper.)

圖3 分解降噪實驗
Fig.3 Experiment de-noising bydecomposition

((a)為原始數(shù)據(jù)；(b)為本文方法的分解降噪結(jié)果。(a)Original data;(b)The result of proposed method.)

圖4是IPIX實測數(shù)據(jù)的降噪實驗結(jié)果，使用了IPIX數(shù)據(jù)文件31#HV極化方式的1 500個樣點，主要用以實驗本文方法的有效性。由圖4可見，本文的分解降噪方法的確較為有效的抑制了噪聲干擾。以此為基礎(chǔ)，可進行降噪數(shù)據(jù)的分析，建模等研究。由于時間限制，這些研究工作有待后續(xù)進行。

3 結(jié)語

信號的統(tǒng)計特性是處理方法選擇的重要依據(jù)。高斯噪聲的相位服從均勻分布，意味著其具有最強的隨機性。因此，本文利用細致的TV方法對矢量序列的相位進行降噪；又根據(jù)矢量的矢徑序列一般服從非均勻分布特點，利用適用于非線性非平穩(wěn)序列的EMD方法進行矢徑分解；據(jù)此形成一種新的矢量序列分解、降噪方法。這種針對不同統(tǒng)計特性“分而治之”的信號處理方法，希望能夠為讀者提供一定的參考。以此為基礎(chǔ)，后續(xù)將進行降噪之后的數(shù)據(jù)分析，建模等研究。