基于CZT的水中目標輻射噪聲的線譜估計方法

李鵬+賈則+張臣+張玉濤

摘要： 基于水下運動目標輻射噪聲線譜數學模型，通過分析水下目標輻射噪聲特征等信息實現。基于經典譜與現代譜估計的聯合估計法提取水下目標輻射噪聲線譜特征，首先估計采集數據的功率譜，得到估算功率譜圖；分析功率譜圖，確定線譜存在的大致區域，再利用線性調頻z變換進行頻譜細化，精確估計線譜位置。

Abstract： Based on the line-spectrum mathematical model of radiation noise pattern of underwater moving target， we can analyse the characteristic of underwater target radiated noise information such as the implementation. Joint estimation method which is based on the classical and modern spectrum estimation can collect the line-spectrum characteristics of underwater target radiation noise. First the power spectrum to estimate the sampling data， and then estimation of power spectrum diagram is obtained. Analysis of power spectrum diagram， we can determine the line-spectrum which is roughly the area， then we can use the chirp z transform frequency spectrum refinement， whicn can accurately estimate the position of line-spectrum.

關鍵詞： 水下運動目標；輻射噪聲；噪聲線譜；CZT

Key words： underwater moving target；radiation noise；line-spectrum of noise；CZT

中圖分類號：TN911.72 文獻標識碼：A 文章編號：1006-4311（2014）15-0027-03

0 引言

線譜是水中目標輻射噪聲功率譜中重要的特征信息，具有強度高、穩定性好以及傳播距離遠等特點。基于線譜分析了解聲源的相關信息，以便進一步識別目標。因此，線譜的檢測在艦艇聲納、魚雷自導以及水中武器引信信號處理中具有舉足輕重的地位。

線譜檢測通常使用周期圖法進行譜估計，分析頻率范圍一般在0～1kHz。對檢測到的線狀譜，要精確測定其頻率并跟蹤測量其變化，提高頻譜的分辨力尤為重要。但在實際應用中，頻譜分辨率與檢測系統處理計算量是一對矛盾，如果使頻譜的采樣集中在線譜附近的頻帶內，就可以在檢測系統處理計算量允許的范圍內盡可能地提高局部頻段內的分辨能力，從而實現對線譜的精確跟蹤測量。鑒于此，筆者分析水中目標輻射噪聲的情況后，提出應用現代譜估計與線性調頻z變換對水下目標輻射噪聲特征線譜進行聯合估測的方法。

1 輻射噪聲線譜數學模型

對于連續譜與線譜所形成的輻射噪聲譜，機械振動及螺旋槳轉動所產生的噪聲均為周期性噪聲源。假設連續譜用寬平穩隨機過程擬合，若考慮調幅特性，可應用局部平穩過程擬合[1，2]。若通過隨機過程G（t）來表示水下運動目標的輻射噪聲，則：

G（t）=L（t）+■X■（t）（1）

式（1）中，L（t）、X■（t）表示寬帶平穩隨機過程和表示初相隨機的周期信號。i=1，…，n說明周期信號有n個。基于此，可通過式（2）計算輻射噪聲功率譜：

G（f）=■■EGk，N（f）■（2）

式（2）中，N為做傅里葉變換時截取的每一段信號長度，E[ ]代表集合平均，k代表信號段的編號。式（2）基于數學模型明確定義了噪聲譜特點，要求獲取無窮多個趨于無窮的長信號段。但是現實情況下只能做有限長以及有限數的集合平均。因此，須通過式（3）確定單邊功率譜：

■（f）=■G■（f）■（3）

假設噪聲信息中包含若干離散頻率成分，則式（1）可寫為：G（t）=L（t）+■a■e■（4）

式（4）中，an屬確定量，即頻率為fn的分量幅度，由此可得出G（t）的Fourier變換公式：

Gk，N（f）=■G（t）e-j2πftdt=Lk，N（f）+■a■■（5）

將變換結果代入式（2）作平均時，設連續譜成分與線譜成分以及不同頻率的線譜成分互相獨立，并令：

C（f）=■■EL■（f）L■■（f） A■=E[α■■]

注意到■■■，在N→∞的情況下，如果f≠fn，它趨于0；如果f=fn，它則趨于∞，因此具有δ函數的性質。繼而得到：

G（f）=C（f）+■Anδ（f-fn）（6）

表明時間信號中的離散頻率成分在功率譜中產生

線譜。

2 水下目標特征譜線聯合估計分析法

2.1 參數模型法

參數模型法是目前較為先進的譜估計方法[3]。AR模型法可整合為計算AR模型各參數ak的課題。筆者擬用Burg遞推算法在Levinson-Durbin遞歸約束的條件下通過計算得出參數ak。Burg遞推算法無需估計自相關函數，它其實是使前向與后向預測誤差能量之和為最小的AR模型功率譜估計方法，能夠保證滿足穩定性的充分條件：k■<1，（p=1，2，…）。而且Burg遞推算法可以準確估計短序列x（n），應用性較強。筆者所提到的水下目標特征譜線聯合估計分析法亦是通過其應用型特點來對特征線譜的位置進行粗略估計，繼而通過線性調頻z變換進行頻譜細化，精確估計線譜位置。endprint

2.2 線性調頻Z變換（CZT）

信號的傅里葉變換處理，可以理解為Z變換在單位圓上進行等間隔采樣。但在很多場合中，并不是整個單位圓上的頻譜都有意義[4]。對于窄帶信號，通常只須分析信號所在的一段頻域范圍，并在這一頻帶內集中采樣，頻帶之外的區域可不作考慮。在這種情況下，可嘗試CZT法分析

頻譜[5]。

假設x（n）為有限長序列，0？燮n？燮N-1，其Z變換為：

X（z）=■x（n）z-n 0？燮n？燮N-1（7）

沿Z平面上的一段螺線進行M點抽樣，抽樣點如下：

zk=AW-k 0？燮k？燮M-1

A、W為復數，二者的極坐標形式為：

A=A■e■ W=W■e■

則zk=A■e■W■e■

式中，A■、W■均為實數；在k=0的情況下，z0=A■e■。由此可見，譜分析起始點z0須根據A0而定，而分析路徑的盤旋趨勢則取決于W0；兩個相鄰分析點之間的夾角是φ0。在A0=r的條件下，W0=1，則須在半徑為r的圓上布設zk的抽樣點；在A0=r=W0=1的情況下，須在單位圓上布設zk的抽樣點，此時相當于DFT變換；當W0<1時，隨著k值增大，zk抽樣點以φ0為步長向外盤旋；W0>1時向內盤旋[6]。通過改變φ0即可獲得不同的頻率分辨率。

將X（zk）改寫為：

X（zk）=■x（n）（AW-k）-n=■x（n）A-nWnk 0？燮k？燮M-1（8）

令nk=■n■+k■-（k-n）■則：

X（zk）=■x（n）A-nW■

=W■■x（n）A■W■W■ 0？燮k？燮M-1（9）

令y（n）=x（n）A■W■，h（n）=W■

則X（zk）改寫為：

X（zk）=W■■y（n）h（k-n） 0？燮k？燮M-1（10）

因此CZT可以利用圖2所示算法框圖實現。

2.3 聯合分析法

筆者基于上述兩種譜估計方法研究出了一種通過聯合分析方法進行水下目標輻射噪聲線譜估計的新方法。具體來講，這種方法是基于AR模型法對x（n）序列估計時不認為N以外的數據全為0，因而估計得到的功率譜中很少有相鄰譜峰的相對值比3dB低的情況，由此便可以得到真實線譜的位置區域，即確定布有特征線譜的區域，繼而通過線性調頻z變換進一步細化頻譜，對線譜位置進行準確預測。

①通過Burg遞推算法估計所采集的水中目標輻射噪聲數據的功率譜，由此得到噪聲線譜的Burg估計算法功率譜圖；

②對①得到的譜估計圖進行分析，根據譜估計圖大致確定線譜的位置；

③基于噪聲線譜大致位置，使用CZT算法對該區域進行頻譜細化，從而對①中的真實線譜區域進行精確線譜估計。

3 仿真研究

基于水下目標輻射噪聲功率譜模型的分析得到輻射噪聲譜的線譜，再用周期信號作模型擬合。假設運動目標相對于水聽器作勻速直線運動，則相對運動速度、水中聲速、水聽器與運動目標航向間的垂直距離以及運動目標和水聽器的連線與其航向的夾角分別用v、c、d和θ表示，t時刻水聽器所接收的、受環境噪聲污染的運動目標輻射噪聲為：

x（t）=■si（t）+n（t）=■Aicos[2π（f0 i+Δfi）（t-ti）+φ0 i]+n（t）（11）

式（11）中：n——不同頻率正弦信號總數；Ai——正弦信號幅度；φ0 i——隨機初相位；f0i——輻射噪聲中的線譜頻率（Hz），則接收信號的多普勒頻率（Hz）可通過Δfi=±■計算得出。在運動目標遠離水聽器接近的情況下取“-”，接近時取“+”，水聽器接收信號的時延（s）通過ti=d/c×sinθ得出，n（t）為噪聲。

假設待分析的信號由水下運動目標的兩個相鄰頻率的正弦信號加噪聲組成，取相鄰頻率分別為f1=150Hz，f2=150.5Hz，θ=40°，c=1500m/s，d=500m，v=40kn，背景噪聲n（t）為高斯白噪聲，其輸入信噪比SNR=5dB。由此對本文的算法進行仿真。

基于AR模型法得到圖3（a）所示的功率譜。根據這一功率譜能夠獲知150Hz周圍的確分布著特征線譜，但不能分辨出150Hz和155Hz兩條譜線。之后再使用CZT算法對150Hz附近的區域進行頻譜細化，這里選取的細化區域為140～160Hz，細化后的頻譜如圖3（b）所示，從圖3（b）中可以精確分辨出150Hz和151Hz兩條譜線。

4 結論

采用周期圖法進行線譜估計，分辨率一般比較低。現代譜估計與信號變換的聯合估計法的應用使得輻射噪聲線譜估計的分辨率大幅提高，有效克服了這一難題。而且，聯合估計法的應用能夠在低信噪比條件下對特征線譜的位置進行準確探測，同時能夠提取精確的特征線譜，工作效率大大提升。

參考文獻：

[1]侯自強.聲納信號處理[M].北京：海洋出版社，1986.

[2]吳國清，李靖，陳耀明等.艦船噪聲識別（I）-總體框架、線譜分析和提取[J].聲學學報，1998，23（5）：394-400.

[3]胡廣書.數字信號處理[M].北京：清華大學出版社，2003.

[4]He Bing， Cabestaing F Postaire， et al. Narrow-band frequency analysis for laser-based lass thickness measurement[J].Instrumentation，2005，54（1）：222-227.

[5]李和明，康偉，顏湘武，張麗霞.一種基于CZT的閃變測量方法[J].電工技術學報，2009，24（3）：209-215.

[6]Rabiner L R，Schafer R W.The Chirp Z transform algorithm[J].IEEE Transactions on Audio and Electroacoustics，1969，17（2）：86-92.endprint

2.2 線性調頻Z變換（CZT）

信號的傅里葉變換處理，可以理解為Z變換在單位圓上進行等間隔采樣。但在很多場合中，并不是整個單位圓上的頻譜都有意義[4]。對于窄帶信號，通常只須分析信號所在的一段頻域范圍，并在這一頻帶內集中采樣，頻帶之外的區域可不作考慮。在這種情況下，可嘗試CZT法分析

頻譜[5]。

假設x（n）為有限長序列，0？燮n？燮N-1，其Z變換為：

X（z）=■x（n）z-n 0？燮n？燮N-1（7）

沿Z平面上的一段螺線進行M點抽樣，抽樣點如下：

zk=AW-k 0？燮k？燮M-1

A、W為復數，二者的極坐標形式為：

A=A■e■ W=W■e■

則zk=A■e■W■e■

式中，A■、W■均為實數；在k=0的情況下，z0=A■e■。由此可見，譜分析起始點z0須根據A0而定，而分析路徑的盤旋趨勢則取決于W0；兩個相鄰分析點之間的夾角是φ0。在A0=r的條件下，W0=1，則須在半徑為r的圓上布設zk的抽樣點；在A0=r=W0=1的情況下，須在單位圓上布設zk的抽樣點，此時相當于DFT變換；當W0<1時，隨著k值增大，zk抽樣點以φ0為步長向外盤旋；W0>1時向內盤旋[6]。通過改變φ0即可獲得不同的頻率分辨率。

將X（zk）改寫為：

X（zk）=■x（n）（AW-k）-n=■x（n）A-nWnk 0？燮k？燮M-1（8）

令nk=■n■+k■-（k-n）■則：

X（zk）=■x（n）A-nW■

=W■■x（n）A■W■W■ 0？燮k？燮M-1（9）

令y（n）=x（n）A■W■，h（n）=W■

則X（zk）改寫為：

X（zk）=W■■y（n）h（k-n） 0？燮k？燮M-1（10）

因此CZT可以利用圖2所示算法框圖實現。

2.3 聯合分析法

筆者基于上述兩種譜估計方法研究出了一種通過聯合分析方法進行水下目標輻射噪聲線譜估計的新方法。具體來講，這種方法是基于AR模型法對x（n）序列估計時不認為N以外的數據全為0，因而估計得到的功率譜中很少有相鄰譜峰的相對值比3dB低的情況，由此便可以得到真實線譜的位置區域，即確定布有特征線譜的區域，繼而通過線性調頻z變換進一步細化頻譜，對線譜位置進行準確預測。

①通過Burg遞推算法估計所采集的水中目標輻射噪聲數據的功率譜，由此得到噪聲線譜的Burg估計算法功率譜圖；

②對①得到的譜估計圖進行分析，根據譜估計圖大致確定線譜的位置；

③基于噪聲線譜大致位置，使用CZT算法對該區域進行頻譜細化，從而對①中的真實線譜區域進行精確線譜估計。

3 仿真研究

基于水下目標輻射噪聲功率譜模型的分析得到輻射噪聲譜的線譜，再用周期信號作模型擬合。假設運動目標相對于水聽器作勻速直線運動，則相對運動速度、水中聲速、水聽器與運動目標航向間的垂直距離以及運動目標和水聽器的連線與其航向的夾角分別用v、c、d和θ表示，t時刻水聽器所接收的、受環境噪聲污染的運動目標輻射噪聲為：

x（t）=■si（t）+n（t）=■Aicos[2π（f0 i+Δfi）（t-ti）+φ0 i]+n（t）（11）

式（11）中：n——不同頻率正弦信號總數；Ai——正弦信號幅度；φ0 i——隨機初相位；f0i——輻射噪聲中的線譜頻率（Hz），則接收信號的多普勒頻率（Hz）可通過Δfi=±■計算得出。在運動目標遠離水聽器接近的情況下取“-”，接近時取“+”，水聽器接收信號的時延（s）通過ti=d/c×sinθ得出，n（t）為噪聲。

假設待分析的信號由水下運動目標的兩個相鄰頻率的正弦信號加噪聲組成，取相鄰頻率分別為f1=150Hz，f2=150.5Hz，θ=40°，c=1500m/s，d=500m，v=40kn，背景噪聲n（t）為高斯白噪聲，其輸入信噪比SNR=5dB。由此對本文的算法進行仿真。

基于AR模型法得到圖3（a）所示的功率譜。根據這一功率譜能夠獲知150Hz周圍的確分布著特征線譜，但不能分辨出150Hz和155Hz兩條譜線。之后再使用CZT算法對150Hz附近的區域進行頻譜細化，這里選取的細化區域為140～160Hz，細化后的頻譜如圖3（b）所示，從圖3（b）中可以精確分辨出150Hz和151Hz兩條譜線。

4 結論

采用周期圖法進行線譜估計，分辨率一般比較低。現代譜估計與信號變換的聯合估計法的應用使得輻射噪聲線譜估計的分辨率大幅提高，有效克服了這一難題。而且，聯合估計法的應用能夠在低信噪比條件下對特征線譜的位置進行準確探測，同時能夠提取精確的特征線譜，工作效率大大提升。

參考文獻：

[1]侯自強.聲納信號處理[M].北京：海洋出版社，1986.

[2]吳國清，李靖，陳耀明等.艦船噪聲識別（I）-總體框架、線譜分析和提取[J].聲學學報，1998，23（5）：394-400.

[3]胡廣書.數字信號處理[M].北京：清華大學出版社，2003.

[4]He Bing， Cabestaing F Postaire， et al. Narrow-band frequency analysis for laser-based lass thickness measurement[J].Instrumentation，2005，54（1）：222-227.

[5]李和明，康偉，顏湘武，張麗霞.一種基于CZT的閃變測量方法[J].電工技術學報，2009，24（3）：209-215.

[6]Rabiner L R，Schafer R W.The Chirp Z transform algorithm[J].IEEE Transactions on Audio and Electroacoustics，1969，17（2）：86-92.endprint

2.2 線性調頻Z變換（CZT）

信號的傅里葉變換處理，可以理解為Z變換在單位圓上進行等間隔采樣。但在很多場合中，并不是整個單位圓上的頻譜都有意義[4]。對于窄帶信號，通常只須分析信號所在的一段頻域范圍，并在這一頻帶內集中采樣，頻帶之外的區域可不作考慮。在這種情況下，可嘗試CZT法分析

頻譜[5]。

假設x（n）為有限長序列，0？燮n？燮N-1，其Z變換為：

X（z）=■x（n）z-n 0？燮n？燮N-1（7）

沿Z平面上的一段螺線進行M點抽樣，抽樣點如下：

zk=AW-k 0？燮k？燮M-1

A、W為復數，二者的極坐標形式為：

A=A■e■ W=W■e■

則zk=A■e■W■e■

式中，A■、W■均為實數；在k=0的情況下，z0=A■e■。由此可見，譜分析起始點z0須根據A0而定，而分析路徑的盤旋趨勢則取決于W0；兩個相鄰分析點之間的夾角是φ0。在A0=r的條件下，W0=1，則須在半徑為r的圓上布設zk的抽樣點；在A0=r=W0=1的情況下，須在單位圓上布設zk的抽樣點，此時相當于DFT變換；當W0<1時，隨著k值增大，zk抽樣點以φ0為步長向外盤旋；W0>1時向內盤旋[6]。通過改變φ0即可獲得不同的頻率分辨率。

將X（zk）改寫為：

X（zk）=■x（n）（AW-k）-n=■x（n）A-nWnk 0？燮k？燮M-1（8）

令nk=■n■+k■-（k-n）■則：

X（zk）=■x（n）A-nW■

=W■■x（n）A■W■W■ 0？燮k？燮M-1（9）

令y（n）=x（n）A■W■，h（n）=W■

則X（zk）改寫為：

X（zk）=W■■y（n）h（k-n） 0？燮k？燮M-1（10）

因此CZT可以利用圖2所示算法框圖實現。

2.3 聯合分析法

筆者基于上述兩種譜估計方法研究出了一種通過聯合分析方法進行水下目標輻射噪聲線譜估計的新方法。具體來講，這種方法是基于AR模型法對x（n）序列估計時不認為N以外的數據全為0，因而估計得到的功率譜中很少有相鄰譜峰的相對值比3dB低的情況，由此便可以得到真實線譜的位置區域，即確定布有特征線譜的區域，繼而通過線性調頻z變換進一步細化頻譜，對線譜位置進行準確預測。

①通過Burg遞推算法估計所采集的水中目標輻射噪聲數據的功率譜，由此得到噪聲線譜的Burg估計算法功率譜圖；

②對①得到的譜估計圖進行分析，根據譜估計圖大致確定線譜的位置；

③基于噪聲線譜大致位置，使用CZT算法對該區域進行頻譜細化，從而對①中的真實線譜區域進行精確線譜估計。

3 仿真研究

基于水下目標輻射噪聲功率譜模型的分析得到輻射噪聲譜的線譜，再用周期信號作模型擬合。假設運動目標相對于水聽器作勻速直線運動，則相對運動速度、水中聲速、水聽器與運動目標航向間的垂直距離以及運動目標和水聽器的連線與其航向的夾角分別用v、c、d和θ表示，t時刻水聽器所接收的、受環境噪聲污染的運動目標輻射噪聲為：

x（t）=■si（t）+n（t）=■Aicos[2π（f0 i+Δfi）（t-ti）+φ0 i]+n（t）（11）

式（11）中：n——不同頻率正弦信號總數；Ai——正弦信號幅度；φ0 i——隨機初相位；f0i——輻射噪聲中的線譜頻率（Hz），則接收信號的多普勒頻率（Hz）可通過Δfi=±■計算得出。在運動目標遠離水聽器接近的情況下取“-”，接近時取“+”，水聽器接收信號的時延（s）通過ti=d/c×sinθ得出，n（t）為噪聲。

假設待分析的信號由水下運動目標的兩個相鄰頻率的正弦信號加噪聲組成，取相鄰頻率分別為f1=150Hz，f2=150.5Hz，θ=40°，c=1500m/s，d=500m，v=40kn，背景噪聲n（t）為高斯白噪聲，其輸入信噪比SNR=5dB。由此對本文的算法進行仿真。

基于AR模型法得到圖3（a）所示的功率譜。根據這一功率譜能夠獲知150Hz周圍的確分布著特征線譜，但不能分辨出150Hz和155Hz兩條譜線。之后再使用CZT算法對150Hz附近的區域進行頻譜細化，這里選取的細化區域為140～160Hz，細化后的頻譜如圖3（b）所示，從圖3（b）中可以精確分辨出150Hz和151Hz兩條譜線。

4 結論

采用周期圖法進行線譜估計，分辨率一般比較低。現代譜估計與信號變換的聯合估計法的應用使得輻射噪聲線譜估計的分辨率大幅提高，有效克服了這一難題。而且，聯合估計法的應用能夠在低信噪比條件下對特征線譜的位置進行準確探測，同時能夠提取精確的特征線譜，工作效率大大提升。

參考文獻：

[1]侯自強.聲納信號處理[M].北京：海洋出版社，1986.

[2]吳國清，李靖，陳耀明等.艦船噪聲識別（I）-總體框架、線譜分析和提取[J].聲學學報，1998，23（5）：394-400.

[3]胡廣書.數字信號處理[M].北京：清華大學出版社，2003.

[4]He Bing， Cabestaing F Postaire， et al. Narrow-band frequency analysis for laser-based lass thickness measurement[J].Instrumentation，2005，54（1）：222-227.

[5]李和明，康偉，顏湘武，張麗霞.一種基于CZT的閃變測量方法[J].電工技術學報，2009，24（3）：209-215.

[6]Rabiner L R，Schafer R W.The Chirp Z transform algorithm[J].IEEE Transactions on Audio and Electroacoustics，1969，17（2）：86-92.endprint