李 欽，劉利民，鞠 峰，韓壯志，何 強

(解放軍軍械工程學院，石家莊 050003)

0 引言

彈丸出膛時刻在火炮動力學分析、實驗測試等領域有著重要作用。延遲測量方法是雷達進行初速測量的常用方法，采用這種方法必須精確給出出膛時刻并作為初速外推的時間零點，可見出膛時刻的精度對初速測量的精度起著至關重要的作用［1－2］。

超高射速火炮彈丸射速快，連發(fā)時炮口火焰噪聲強，傳統(tǒng)的出膛時刻［3－4］獲取方法已不能滿足超高射速彈丸出膛時刻獲取的要求。由于超高射速武器彈丸出膛時的回波較弱，有效回波淹沒在噪聲中，采用雷達測速只能獲取出膛時刻的粗略估計，無法進行精確定位，嚴重影響了初速測量的精度。

針對上述問題，文中對傳統(tǒng)的雷達測速方法進行改進，在時頻分析的基礎上增加降噪處理，用圖像邊緣檢測的方法實現對回波信號突變點的精確定位，獲得超高射速彈丸的出膛時刻。

1 算法原理

1.1 短時傅里葉變換

短時傅里葉變換是處理非平穩(wěn)信號的經典方法，通過時間域上的加窗處理，將非平穩(wěn)信號分解成一系列短時平穩(wěn)信號，并分別對每一個短時平穩(wěn)信號進行處理，最后通過時間參數t在整個時間域的平移完成對非平穩(wěn)信號的處理［5］。工程實際通常需要將STFT(t，f)離散化，其短時譜定義為:

式中:g(μ)為歸一化窗函數;λ是分幀的序列號;k為頻率點號，k=0，1，…，L－1;L為窗函數寬度;R為幀移長度。S(λ，k)是信號在時間窗窄區(qū)間內的頻譜，又稱為短時譜。得到正確有效的短時譜，是進行時頻分析和基于時頻分析降噪的前提。

1.2 最小統(tǒng)計噪聲譜估計

最小統(tǒng)計噪聲譜估計首先在時間窗口內搜索子帶功率的最小值，然后把經過誤差補償的最小值作為該頻帶的噪聲譜估計，最后利用譜減法得到較為純凈的信號功率譜。

假設 Y(λ，k)為雷達回波信號短時譜，S(λ，k)為彈丸回波信號短時譜，N(λ，k)為噪聲短時譜，則有:

在長度為D的滑動窗口內對雷達回波信號功率譜進行搜索，找到窗內每個頻點的功率最小值［6－7］:

隨機變量的最小值小于其均值，所以需要對得到的最小功率值Ymin(λ，k)進行補償，這里引入一個偏差補償因子omin(D)。Yd(λ，k)是噪聲估計功率譜:

式中，偏差補償omin(D)是一常數，omin(D)取值越大，降噪效果越明顯，但同時也帶來了信號的失真。

利用譜減法思想［8］，彈丸回波信號的功率譜S(λ，k)可以通過雷達回波信號的功率譜Y(λ，k)與估計的噪聲功率譜Yd(λ，k)相減得到:

基于最小統(tǒng)計噪聲譜估計的譜減法，能對噪聲的估計值及時更新，殘留噪聲大大減少，降噪效果比較明顯。

1.3 圖像邊緣檢測

圖像邊緣是指像素灰度發(fā)生階躍變化或屋頂狀變化的像素點的集合［9－10］。圖像邊緣檢測是對灰度變化的度量與定位，而灰度變化可以通過導數來度量，因此邊緣檢測的一個基本思想就是求解一階導數的局部極大值或者二階導數的過零點。

利用梯度最大值提取邊緣點的這種思想產生了許多經典的邊緣檢測方法，如Sobel算法、Log算法、Canny算法等。其中，Sobel算法結合了高斯平滑和微分，相比于其他算法具有邊緣提取效果好，對噪聲具有一定魯棒性，運算簡單，易于硬件實現且實時性好的優(yōu)點［11］。因此，文中選取了Sobel算法對圖像進行邊緣檢測。

1.4 Sobel算子和二值化

Sobel算子是基于一階導數的邊緣檢測算子，求出的邊緣點存在于圖像梯度最大值處。梯度是函數變化的一種度量，梯度幅度值和方向定義為:

上式給出了灰度在x和y方向上的變化率，若計算灰度時均采用上式，運算量較大。通常使用小型差分模板利用卷積來近似計算。Sobel算子中用到了兩個卷積模板，前一個對水平邊緣影響較大，稱為水平算子;后一個對垂直邊緣影響較大，稱為垂直算子。

Sobel算子是全方向微分算子，對檢測點的上下左右進行加權，但鄰點對當前點產生的影響不是等價的，距離不同的點，權值不同，對算子結果的影響也不同［10］。Sobel算子先用兩個卷積模板對圖像進行卷積，然后根據梯度公式得到各個像素灰度梯度值的平方并與閾值進行比較。在閾值不為空的情況下，邊緣存在于梯度幅值平方大于閾值平方的點上。

Sobel算法采用單閾值判決，閾值的選取至關重要?？紤]到實際情況采用平均值加權的方法選取閾值。對圖片內的像素值求和取平均，再乘以加權因子作為閾值。

所謂的圖像二值化就是將多值的灰度圖像轉化為“0”和“1”表示的二值圖像。在文中指將邊緣檢測后的圖形內部和邊緣的像素賦值為“1”，其他部分像素賦值為“0”，方便后續(xù)出膛時刻的檢測和判定。

2 算法流程及結果分析

2.1 算法流程

算法的實現主要包括以下步驟:

1)通過短時傅里葉變換獲取時頻二維矩陣;

2)最小統(tǒng)計噪聲譜估計法進行降噪處理;

3)將回波對數能量以二維矩陣形式存儲，并保存為一張圖片;

4)對圖片用Sobel算子進行邊緣檢測，提取邊緣信息并二值化處理;

5)確定每發(fā)彈丸回波的起始點。

首先對接收到的回波數據進行短時傅里葉變換，把一維時域信號轉化為二維時頻域信號。由于回波信號有較強的火焰噪聲，變換后得到的頻譜圖像不能清晰的反映彈丸回波的起始點，必須進行降噪處理，而“最小統(tǒng)計噪聲譜估計”方法，可以在很大程度上抑制火焰噪聲。

將降噪后的回波能量保存為一張圖片的形式，那么彈丸的回波能量分布就反映為圖片中圖形的分布情況。然后用1.4中提到的Sobel算法進行邊緣檢測，提取圖片的邊緣信息，并進行二值化處理，方便后續(xù)的出膛時刻及彈丸個數確定。

在處理圖片中彈丸信息時，每個圖形的列數代表每發(fā)彈丸的持續(xù)時間，最小的列號即彈丸出膛時刻，不同的彈丸持續(xù)時間不同;行數代表彈丸回波信號能量的幅度，最小的行號即彈丸信號峰值。兩個行數閾值可以確定圖形是否有效:低閾值與最小行號比較，確定信號是否有效;高閾值用來消除圖形斷裂帶來的影響。

2.2 結果分析

本次仿真采用實測的步槍信號模擬高射速火炮信號，系統(tǒng)的采樣周期T=1 μs，系統(tǒng)的時間精度要求T0=0.1 ms。圖1為步槍在不均勻射擊的情況下6連發(fā)彈丸的實測雷達回波時域信號。橫軸表示時間，縱軸表示幅度。從信號的時域圖中可以明顯看出信號有6個峰值，表明有6發(fā)槍彈，但是由于存在較大的噪聲，無法精確地在時域上直接檢測彈丸的出膛時刻。

圖1 回波信號時域圖

圖2是短時傅里葉變換后的時頻圖，橫軸代表時間(幀數)，縱軸代表頻率，明亮程度代表幅度的大小，越亮表示幅度值越大。這里采用blackman-harris窗作為窗函數，窗寬256，幀移100。由于FFT變換具有對稱性，所以畫出前128個頻點。6根明亮譜線開始的時刻則對應著每發(fā)彈丸的出膛時刻。短時譜圖由于強烈火焰噪聲的影響，無法準確檢測出彈丸的出膛時刻，需要進行降噪處理。

圖2 回波信號時頻圖

圖3是降噪之后的時頻圖，橫軸代表時間(幀數)，縱軸代表頻率，由于能量大都集中在低頻部分，所以只取了前32個頻點的值。采用最小統(tǒng)計噪聲譜估計法降噪時，搜索窗的長度D必須包含一個完整的彈丸回波信號，考慮到精度要求又不宜太小，設置D=20;根據經驗值和實際效果，偏差補償omin(D)=4。從圖中可以看出降噪效果比較明顯，濾除了火焰噪聲，彈丸回波圖形能量分布清晰可辨，為起始點檢測奠定了基礎。

圖3 降噪后回波信號時頻圖

圖4 邊緣檢測二值化的圖形

圖4是經過Sobel邊緣檢測并二值化處理的結果圖。由于圖片的背景反差比較大，通過實驗發(fā)現閾值加權系數取2能得到比較理想的結果。圖像黑色的區(qū)域像素值為0，白色區(qū)域像素值為1。每一個白色區(qū)域表示1發(fā)彈丸，從圖中可以清晰地看出有6發(fā)彈丸，便于直接進行邊緣起始點的檢測。

得到圖4圖形之后，對每一個白色圖形區(qū)域進行峰值檢測和邊緣起始點的提取。在處理圖片中彈丸信息時，每個圖形的列數代表每發(fā)彈丸的持續(xù)時間，根據精度要求持續(xù)時間不得小于0.5 ms，即列數不得少于5列，得到的最小的列號即該發(fā)彈丸的出膛時刻。行數代表彈丸回波信號能量的幅度，最小的行號即彈丸信號峰值，最小行號用來判斷該圖形是否為有效彈丸，滿足兩個閾值的為有效彈丸。低閾值為最大行數減去70，高閾值為最大行數減去10。

這組數據得到的結果為:6發(fā)彈丸，每一發(fā)彈丸的出膛時刻:15.0 ms、18.9 ms、54.2 ms、114.9 ms、129.9 ms、150.8 ms。從最后得到的出膛時刻可以清晰的檢測出發(fā)射間隔為3.9 ms的彈丸，即射速大約15 000發(fā)/min的高射速彈丸。這說明本算法能夠滿足系統(tǒng)對于精度的要求，同時也驗證了本算法精確獲取高射速火炮彈丸出膛時刻的可行性和有效性，為測速雷達解算彈丸初速提供了時間零點。

3 結論

文中提出的圖像邊緣檢測的處理方法，能夠直接在時頻域對彈丸回波進行處理，并精確定位和獲取彈丸的出膛時刻，最后通過實測數據進行了驗證。實測數據表明該算法實現簡單，可以精確的獲取彈丸的出膛時刻，為今后在獲取彈丸出膛時刻的方法上提供了參考。

［1］胡江，黃景徳，解維河.基于測速雷達的艦炮初速測量技術研究［J］.艦船電子工程，2011，31(6):94-96.

［2］馬玲，蔡征宇，程風雷，等.毫米波測速雷達的測速原理［J］.彈道學報，2003，15(4):87-91.

［3］王寶元，鈔紅曉，邵小軍，等.彈丸出炮口時間測試方法研究［J］.兵工學報，2012，33(6):736-740.

［4］劉海林.獲取測速雷達時間零點的幾種方法［J］.無線電工程，2008，38(3):61-64.

［5］張賢達，保錚.非平穩(wěn)信號分析與處理［M］.北京，國防工業(yè)出版社，1998.

［6］牛銅，張連海，屈丹.基于加權最小統(tǒng)計的噪聲普估計改進算法［J］.電子與信息學報，2009，31(5):1166-1169.

［7］Martin R.Noise power spectral density estimation based on optimal smoothing and minimum statistics［J］.IEEE Trans.Speech Audio Process，2001，9(5):504-512.

［8］李銀國，薛雯，徐洋.基于噪聲短時譜動態(tài)估計的語音增強譜減算法［J］.重慶郵電大學學報，2010，22(2):127-130.

［9］崔建軍，詹世富，鄭雄偉，等.一種改進的邊緣檢測算法［J］.測繪科學，2009，34(4):55-56.

［10］李婭婭，李志潔，鄭海旭，等.圖像邊緣檢測算法的比較與實現［J］.計算機工程與設計，2010，31(9):1971-1975.

［11］佟文君，張軍.基于邊緣檢測與傳統(tǒng)差分法的目標檢測方法［J］.天津職業(yè)技術師范大學學報，2011，21(2):37-40.