劉 丹，張丕狀，馬春燕

（中北大學(xué) 信息探測(cè)與處理技術(shù)研究所，山西 太原 030051）

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

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

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

1 基于最小二乘法的彈丸回波信號(hào)估計(jì)算法

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

1.1 最小二乘法簡(jiǎn)介

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

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

使得

成立，則稱曲線y*（x）為在曲線族中按最小二乘確定的對(duì)于數(shù)據(jù)（xi，yi）的擬合曲線。 由線性無關(guān)向量組 φj=（φj（x0），φj（x1），L，φj（xm））T（j=0～n）做基底構(gòu)成一個(gè) Rm+1的一個(gè)子空間，記 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 采用最小二乘法擬合將信號(hào)去噪

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

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

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

1.3 擬合方法的改進(jìn)

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

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

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

其中，

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

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

2 基于STFT時(shí)頻分析法分析彈丸回波信號(hào)

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

2.1 STFT時(shí)頻分析法

Potter等在1947年首次提出了短時(shí)Fourier變換[6]。基本思想是，用一個(gè)在時(shí)間軸上可移動(dòng)的、時(shí)間跨度很小的分析窗函數(shù)h（t）截取信號(hào)，并認(rèn)為窗內(nèi)的信號(hào)是準(zhǔn)平穩(wěn)的，進(jìn)行Fourier變換提取頻率分量的時(shí)域局部化信息，然后移動(dòng)窗函數(shù)，重復(fù)上述過程，則得到短時(shí)Fourier變換，也稱為窗口Fourier變換，其表達(dá)式如下：

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

2.2 工程應(yīng)用中的改進(jìn)

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

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

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

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

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

因此，最佳窗長(zhǎng)由下式給出：

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

信號(hào)相應(yīng)的最佳短時(shí)傅里葉變換就可以由下式得到：

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

3 驗(yàn)證分析

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

3.1 對(duì)最小二乘法數(shù)據(jù)擬合的驗(yàn)證

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

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

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

圖1 雷達(dá)采集到的彈丸回波信號(hào)Fig.1 The projectile radar singal collected by radar

圖2 擬合后信號(hào)的對(duì)比圖Fig.2 The signal contrast figure after ftting

圖3 一般擬合方法處理后的信號(hào)對(duì)比圖Fig.3 The signal contrast figure after the general fttingmethod

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

3.2 采用窗長(zhǎng)自適應(yīng)的STFT法對(duì)信號(hào)時(shí)頻分析的對(duì)比

圖5是對(duì)原始信號(hào)進(jìn)行短時(shí)傅里葉變換分析的時(shí)頻圖，圖6是對(duì)回波信號(hào)進(jìn)行擬合處理后，并采取自適應(yīng)短時(shí)傅里葉變換分析的時(shí)頻圖。

圖4 改進(jìn)擬合方法處理后的信號(hào)對(duì)比圖Fig.4 The signal contrast figure after the improved fttingmethod

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

圖6 自適應(yīng)STFT時(shí)頻分析Fig.6 Time-frequency analysis of the adaptive STFT

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

4 結(jié) 論

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

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

[2]王濤，孟凡順，李洋森，等.不同時(shí)頻分析方法的精度比較及應(yīng)用[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]曹健，林濤，徐遐齡，等.基于最小二乘法和時(shí)頻原子變換的諧波/間諧波測(cè)量算法[J].電工技術(shù)學(xué)報(bào)，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].振動(dòng)與沖擊，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]史林，鞠峰，胡文華，等.基于短時(shí)傅里葉變換的高射速火炮彈丸出膛時(shí)刻測(cè)試方法[J].火炮發(fā)射與控制學(xué)報(bào)，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)信號(hào)特征提取方法及其應(yīng)用[M].北京：科學(xué)出版社，2013.

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

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