基于最小二乘法和STFT的彈丸回波信號估計分析

劉 丹，張丕狀，馬春燕

（中北大學 信息探測與處理技術研究所，山西 太原 030051）

炮管膛內(nèi)如運動行程、速度、加速度與時間的關系等參數(shù)的測量對于鑒定武器彈藥特性、槍炮管身特性，評價武器系統(tǒng)壽命、安全性和可靠性具有重要意義[1]。彈丸膛內(nèi)運動是一個短時變速過程，穿越雷達波束的回波信號是頻率隨時間快速連續(xù)變化的非平穩(wěn)信號，頻譜表現(xiàn)為高度密集分布，多普勒頻率動態(tài)范圍也很大，既包含了彈丸在膛內(nèi)的運動參數(shù)也包含了槍炮管身的運動參數(shù)。所以彈丸回波信號的時頻分析就成為研究彈丸在膛內(nèi)運動特征的重要對象。

目前時頻分析方法已成為分析彈丸回波信號廣泛使用的技術，主要包括瞬時屬性提取、信號頻譜分解、信號能量衰減補償、信號識別等應用技術[2]。文中對彈丸回波信號分析采用了短時傅里葉變換。由于彈丸在運行過程中，管身的震動等外界因素使彈丸回波信號產(chǎn)生嚴重的噪聲，所以直接使用短時傅里葉變換對彈丸回波信號進行時頻分析特征不明顯，傳統(tǒng)的小波去噪、數(shù)字濾波會使信號處理后丟失有效信息，影響后續(xù)彈丸速度測量的精度。

文中提出了一種方法對彈丸回波信號進行預處理，消除高頻噪聲。首先在時域內(nèi)對彈丸回波信號數(shù)據(jù)采用最小二乘法進行分段擬合處理，并對擬合方式進行改進，使擬合后的數(shù)據(jù)更加逼近彈丸回波信號的實際情況。然后通過短時傅里葉變換對原始信號和處理后的信號進行時頻分析，通過比較，處理后的信號時頻特征更加清晰，驗證了方法的可行性。

1 基于最小二乘法的彈丸回波信號估計算法

彈丸高速發(fā)射時，由于背景噪聲的存在嚴重影響回波信號的特性，頻譜分析時需要對回波信號進行預處理。曲線擬合是圖像分析中非常重要的描述符號，最常用的曲線擬合方法是最小二乘法。

1.1 最小二乘法簡介

最小二乘法（Least Squares，LS）是一種數(shù)學優(yōu)化技術，其基本思想是:給定一組有序數(shù)對，根據(jù)誤差平方和最小化原則，找出這些數(shù)據(jù)的最佳函數(shù)匹配。

最小二乘法的數(shù)學原理為[3]：假使曲線

使得

成立，則稱曲線y*（x）為在曲線族中按最小二乘確定的對于數(shù)據(jù)（xi，yi）的擬合曲線。 由線性無關向量組 φj=（φj（x0），φj（x1），L，φj（xm））T（j=0～n）做基底構(gòu)成一個 Rm+1的一個子空間，記 A=（φ0，φ1，L，φn），y=（y0，y1，L，ym）T。向量組滿足條件（2）的擬合曲線y*（x）存在且唯一，并且從方程

中解出 c*=（c*0，c*1，…，c*n）即可得到擬合曲線（1.1）。

1.2 采用最小二乘法擬合將信號去噪

已經(jīng)分析得知，彈丸回波信號是頻譜高度密集分布的非平穩(wěn)信號，而且彈丸在膛內(nèi)運行過程中，受震動、火焰等背景噪聲的影響，有效回波信號變?nèi)酰€圖像變得粗糙。為了將回波信號中不必要的噪聲信息去除，本文提出了采用最小二乘法對信號數(shù)據(jù)進行擬合處理，逼近有效的回波信號特征。彈丸回波信號的基本模型符合正弦波的信號特征。對采集到的雷達數(shù)據(jù)，文中采取分段擬合，將信號波段用拋物線逼近擬合，按照樣點在均值線兩側(cè)分布均勻和多項式階數(shù)盡可能低的準則[4]，采取二次多項式函數(shù)進行分段擬合。通常認為下式給出的二次多項式模擬程度較好。

其中 u（t）是每段數(shù)據(jù)擬合出的拋物線，ai（i=0，1，2）是擬合多項式的系數(shù)。然后將擬合出的數(shù)據(jù)段相加即可認為逼近彈丸回波信號的有效信號xt。

n為信號被分成的段數(shù)，依據(jù)信號中每段類拋物線點數(shù)和信號總長度計算確定。

1.3 擬合方法的改進

一般情況下，分段擬合的方法是將信號數(shù)據(jù)分成若干段數(shù)，由于分段的函數(shù)模型差異，每段端點均沒有重合，為此將擬合后的數(shù)據(jù)直接相加去逼近所需分析的信號。文中試驗采集到的數(shù)據(jù)，每兩千個點可構(gòu)成一段拋物線，故將信號數(shù)據(jù)每兩千個點分成一段。例如前兩段，x1（1:2000），x2（2001:4000）擬合后數(shù)據(jù)為，u1（1:2000），u2（2001:4000）則 xt=u1+u2。 后面每段的劃分及處理以此類推。

這樣擬合相加出現(xiàn)一個問題，即在分段信號的相加點處出現(xiàn)一個凸起點，特別是影響單調(diào)性的凸起點，對后續(xù)的信號特征提取產(chǎn)生嚴重影響。

改進的辦法是，每段仍按2 000個點數(shù)據(jù)進行擬合，但從下一段開始，起點從上一段擬合好的數(shù)據(jù)后若干個點處算起，仍舉前兩段為例，則x1=x（1:2 000），第二段起點從第一段擬合后的數(shù)據(jù)后 600 個點處開始，x2=u1（1 401:2 000）+x（2 001:3 400），后面每段的劃分以此類推。從第二段起，每段數(shù)據(jù)可用如下表達式表示：

其中，

上式中，i從1開始計起，M是每段數(shù)據(jù)從上一段擬合后的數(shù)據(jù)中取出后面的點數(shù)，ui是前一段擬合后的數(shù)據(jù)，ui是從中取出的后M個點的數(shù)據(jù)，x是原始數(shù)據(jù)，x′是原始數(shù)據(jù)中取出的2 000－M個點的數(shù)據(jù)，xi+1是從第二段起每段要擬合的總數(shù)據(jù)。則擬合后信號總數(shù)據(jù)為下式：

這樣，經(jīng)過改進擬合方法后對數(shù)據(jù)進行處理，可大大緩解擬合端點處的凸起問題，有效避免了后續(xù)信號特征提取時出現(xiàn)的誤差。

2 基于STFT時頻分析法分析彈丸回波信號

短時傅里葉變換是分析非平穩(wěn)信號的有力工具[5]。文中基于短時傅里葉變換法，分析原信號和數(shù)據(jù)擬合后信號的時頻譜分布，驗證信號預處理方法的可行性。

2.1 STFT時頻分析法

Potter等在1947年首次提出了短時Fourier變換[6]?；舅枷胧?，用一個在時間軸上可移動的、時間跨度很小的分析窗函數(shù)h（t）截取信號，并認為窗內(nèi)的信號是準平穩(wěn)的，進行Fourier變換提取頻率分量的時域局部化信息，然后移動窗函數(shù)，重復上述過程，則得到短時Fourier變換，也稱為窗口Fourier變換，其表達式如下：

短時Fourier變換克服了一般Fourier變換譜分析中時間域無限大的缺點，給信號加一個窗，信號與對應于某一時移和頻移的窗函數(shù)的內(nèi)積就能反映信號的局部頻譜特性，整個變換結(jié)果也就能揭示信號頻譜的演化特性。

2.2 工程應用中的改進

對工程實際應用而言，需將短時傅立葉變換離散化。對于離散序列信號x（m），它在時刻的短時數(shù)據(jù)可以定義為xm（n）=x（n）w（m－n），其中 w（n）為窗函數(shù)，窗的長度為 Nw。 一般選中心對稱的滑動窗，能量主要集中在n或原點（n=0）附近。此時的窗長度為奇數(shù)，使得窗的中心點正好對應要分析的時間點。其離散Fourier變換可表示為：

上式即為非平穩(wěn)信號x（m）的離散STFT變換公式，它給出了信號在m=n附近的一段時間內(nèi)的時頻信息[7]。

眾所周知，時頻分辨率只與窗函數(shù)的時頻跨度有關，而窗函數(shù)的時頻跨度取決于窗函數(shù)形式和長度，短時傅里葉變換采用了固定窗函數(shù)，故其只具有單一分辨率。運用時，須近似認為彈丸回波信號在窗口內(nèi)是平穩(wěn)信號疊加噪聲。

已分析知彈丸出膛是一個短時高速過程，回波信號波形變化劇烈并參有強烈的振動噪聲，不同時刻回波信號的時頻分辨率不盡相同。用窗函數(shù)不變的短時傅里葉變換分析回波信號，會產(chǎn)生時間或頻率上的模糊現(xiàn)象。因此文中采用文獻[8]提出的自適應算法，按照信號歸一化后局部時頻能量最大的原則，為信號時域點逐個選擇最佳窗長，即讓窗長自適應化。

假設輸入信號為 x（t），窗函數(shù)為 g（t），計算每一個信號點短時傅立葉變換后的局部能量，它是以窗長T為自變量的函數(shù)[1]。

因此，最佳窗長由下式給出：

其中 en（k）|x（k+m）|。AT為窗函數(shù)能量的倒數(shù)，作為窗長增加的懲罰因子。<…>表示內(nèi)積。

信號相應的最佳短時傅里葉變換就可以由下式得到：

經(jīng)過改進，再分析回波信號的頻譜，分辨率得到有效提升，體現(xiàn)在時頻圖上就是時頻能量聚集性，聚集性越高分辨率越高。

3 驗證分析

圖1是雷達采集到的彈丸回波信號，文中只截取其中一段來說明。橫坐標為數(shù)據(jù)點數(shù)。我們需要的有效信息是圖中類似正弦波的信號，但從中可以看到，信號的包絡上含有很多噪聲，需做的工作即將這些噪聲最大程度的消除并用短時傅里葉變換分析其頻譜。

3.1 對最小二乘法數(shù)據(jù)擬合的驗證

圖2是經(jīng)過最小二乘法數(shù)據(jù)擬合后的信號。

從對比圖上可以看出，經(jīng)過擬合處理后的信號，邊緣數(shù)據(jù)得到有效剔除，信號信息得到改善。

下面驗證分析擬合方法改進前后的信號圖像特征。不妨將信號數(shù)據(jù)所取點數(shù)進一步縮小，以觀察更為明顯。圖3和圖4截取2 000～10 000個數(shù)據(jù)點的信號圖。紅色為擬合后的曲線。圖3采用一般的擬合方法對數(shù)據(jù)進行處理，圖4采用改進擬合方法對數(shù)據(jù)進行處理。

圖1 雷達采集到的彈丸回波信號Fig.1 The projectile radar singal collected by radar

圖2 擬合后信號的對比圖Fig.2 The signal contrast figure after ftting

圖3 一般擬合方法處理后的信號對比圖Fig.3 The signal contrast figure after the general fttingmethod

可以看到，采用一般擬合方法對信號處理后，圖像在每段的擬合端點處有明顯的凸起點，特別影響圖像的單調(diào)性分析。而采用改進的擬合方法處理后，擬合端點處的凸起點得到明顯消除。

3.2 采用窗長自適應的STFT法對信號時頻分析的對比

圖5是對原始信號進行短時傅里葉變換分析的時頻圖，圖6是對回波信號進行擬合處理后，并采取自適應短時傅里葉變換分析的時頻圖。

圖4 改進擬合方法處理后的信號對比圖Fig.4 The signal contrast figure after the improved fttingmethod

圖5 一般STFT時頻分析Fig.5 Time-frequency analysis of general STFT

圖6 自適應STFT時頻分析Fig.6 Time-frequency analysis of the adaptive STFT

可以看出，擬合后的信號通過自適應STFT法分析后，時頻譜更加清晰。進一步驗證了最小二乘法對信號數(shù)據(jù)擬合的可行性。

4 結(jié) 論

由于彈丸回波信號噪聲的存在，直接使用時頻分析法分析頻譜特征不清晰。提出了基于最小二乘法對信號數(shù)據(jù)進行擬合處理，并針對數(shù)據(jù)接口處的凸起問題，改進了擬合的方法，效果比較明顯。然后通過窗長自適應的STFT法對原始信號和擬合處理后的信號做了時頻分析對比，驗證了方法在工程應用上的可行性。

[1]張瑋.數(shù)字下變頻技術在彈丸回波信號處理中的研究與應用[D].西安：西安電子科技大學，2012.

[2]王濤，孟凡順，李洋森，等.不同時頻分析方法的精度比較及應用[J].海洋地質(zhì)前言，2013，29（3）:60-64.WANG Tao，MENG Fan-shun，LIYang-sen，et al.The precision and application of themethod of different time domain analysismethods[J].Marine Geology Preface，2013，29（3）:60-64.

[3]曹健，林濤，徐遐齡，等.基于最小二乘法和時頻原子變換的諧波/間諧波測量算法[J].電工技術學報，2011，26（10）:1-6.CAO Jian，LIN Tao，XU Xia-ling，et al.Based on the least square method and measuring of the harmonic wave in the time-frequency atom transformation between harmonic/algorithm[J].Transactions of China Electrotechnical Society，2011，26（10）:1-6.

[4]胡燦陽，陳清軍.基于EMD和最小二乘法的基線飄移研究[J].振動與沖擊，2010，29（3）:162-167.HUCan-yang，CHEN Qing-jun.The baseline driftstudy based on the baseline drift of the EMD and least-squaremethod[J].Journal of Vibration and Shock，2010，29（3）:162-167.

[5]史林，鞠峰，胡文華，等.基于短時傅里葉變換的高射速火炮彈丸出膛時刻測試方法[J].火炮發(fā)射與控制學報，2014，35（2）:35-38.SHILin，JU Feng，HUWen-hua，et al.Rof artillery projectile chamber testmethod based on short-time fourier transform[J].Journal of Artillery Launch and Control，2014，35（2）:35-38.

[6]范虹.非平穩(wěn)信號特征提取方法及其應用[M].北京：科學出版社，2013.

[7]萬永革.數(shù)字信號處理的MATLAB實現(xiàn)[M].北京：科學出版社，2012.

[8]CzerwinskiRN.Adaptive short-time Fourier analysis[J].IEEE Signal Processing Letters，1997，4（2）:42-45.