姬婷婷， 張 杰， 王國(guó)宇

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

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

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

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

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

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

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

1 矢量時(shí)間序列的分解

1.1 相位序列TV降噪

任意一個(gè)矢量時(shí)間序列可以表示為

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

(1)

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

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

φ=y+n。

(2)

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

(3)

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

(4)

(5)

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

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

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

1.2 修改的EMD矢徑分解

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

EMD分解框架：

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

(3)局部均值曲線：對(duì)每一個(gè)樣點(diǎn)，均值曲線：

(6)

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

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

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

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

1.3 矢量重構(gòu)

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

(7)

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

1.4 瞬時(shí)相位梯度與瞬時(shí)功率譜計(jì)算

(8)

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

(9)

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

2 分解實(shí)驗(yàn)

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

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

實(shí)驗(yàn)用例2，線性頻率調(diào)制+正弦頻率調(diào)制信號(hào)，過(guò)程同實(shí)驗(yàn)用例1。

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

((a)純凈信號(hào)的瞬時(shí)頻率；(b)噪聲污染信號(hào)的瞬時(shí)頻率；(c)本文方法降噪信號(hào)的瞬時(shí)頻率；(d)純凈信號(hào)的瞬時(shí)功率；(e)噪聲污染信號(hào)的瞬時(shí)功率；(f)本文方法降噪信號(hào)的瞬時(shí)功率。(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 分解降噪實(shí)驗(yàn)

Fig.2 Experiment de-noising by decomposition

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

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

((a)純凈信號(hào)的瞬時(shí)頻率；(b)噪聲污染信號(hào)的瞬時(shí)頻率；(c)本文方法降噪信號(hào)的瞬時(shí)頻率；(d)純凈信號(hào)的瞬時(shí)功率；(e)噪聲污染信號(hào)的瞬時(shí)功率；(f)本文方法降噪信號(hào)的瞬時(shí)功率。(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 分解降噪實(shí)驗(yàn)
Fig.3 Experiment de-noising bydecomposition

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

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

3 結(jié)語(yǔ)

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