朱建軍魏玉闊杜偉東李海森*

①(哈爾濱工程大學水聲技術重點實驗室 哈爾濱 150001)

②(哈爾濱工程大學水聲工程學院 哈爾濱 150001)

基于分數(shù)階傅里葉變換的Chirp淺剖精細探測方法

朱建軍①②魏玉闊①②杜偉東②李海森*①②

①(哈爾濱工程大學水聲技術重點實驗室 哈爾濱 150001)

②(哈爾濱工程大學水聲工程學院 哈爾濱 150001)

實現(xiàn)弱回波信號檢測和高信噪比(SNR)淺剖圖像獲取是淺剖精細探測的首要任務。該文在分析分數(shù)階傅里葉變換(FrFT)解卷積原理，推導時間量綱化變換公式的基礎上，提出一種基于FrFT的淺剖精細探測新方法。該方法通過FrFT解卷積實現(xiàn)分數(shù)階傅里葉域(u域)沉積層沖激響應求解，采用u域加窗濾波技術對帶內(nèi)噪聲進行有效抑制，經(jīng)時間量綱化變換實現(xiàn)高信噪比u域沉積層沖激響應包絡信號至時域淺剖包跡的直接變換，得到高質(zhì)量的淺剖圖像。仿真實驗和實測數(shù)據(jù)處理驗證了算法的精細探測能力，算法性能優(yōu)于脈沖壓縮和自回歸(AR)預測濾波方法。

信號處理；淺地層剖面；精細探測；分數(shù)階傅里葉變換；帶內(nèi)噪聲抑制；時間量綱化變換

1 引言

海底淺地層剖面探測(簡稱淺剖探測)的精細程度直接關系著海底科學考察的有效性、沉積層結(jié)構分析的合理性和遙測參數(shù)獲取的準確性[1]。目前，淺剖探測廣泛采用基于脈沖壓縮技術的Chirp聲吶探測方法[2,3]，同時也將解卷積技術用于信號處理[4]。但上述技術方法在進行淺剖探測時存在以下問題：首先，脈沖壓縮子波通常存在較高的旁瓣，當子波旁瓣與弱子波主瓣幅度相當時會產(chǎn)生子波遮蔽效應，造成地層誤判和漏判。雖然采用探測信號加窗法可抑制脈沖壓縮子波的旁瓣[5]，但窗函數(shù)會減弱信號能量、減小信號有效帶寬，降低探測分辨率和探測深度。其次，在不損失分辨率的前提下有效抑制帶內(nèi)噪聲困難。傳統(tǒng)去噪方法通常具有濾波的性質(zhì)，在濾除噪聲的同時會損失信號帶寬，導致分辨率下降。特別是在對寬帶信號處理時，傳統(tǒng)帶通濾波技術難以濾除帶內(nèi)噪聲干擾[6]，而帶內(nèi)噪聲的存在會影響對弱回波信號的檢測，限制系統(tǒng)探測性能和淺剖成圖質(zhì)量的進一步提高。最后，解卷積處理需要在較高信噪比下才具有良好的處理性能[7,8]，低信噪比時實現(xiàn)精細探測困難。

由于Chirp信號在最優(yōu)階分數(shù)階傅里葉u域上具有最佳的能量聚集性[9]，而噪聲仍均勻分布于幾乎整個分數(shù)階傅里葉域[10]?；谶@一優(yōu)越處理性能，本文提出一種基于分數(shù)階傅里葉變換(Fractional Fourier Transform, FrFT)技術的Chirp淺剖精細探測方法，通過FrFT解卷積、u域加窗濾波和時間量綱化變換，在u域?qū)崿F(xiàn)高信噪比淺剖圖像的獲取和淺地層剖面的更精細探測，達到傳統(tǒng)時域或頻域方法均無法獲得的處理效果，仿真實驗和實測數(shù)據(jù)處理驗證了方法的有效性。

2 淺剖回波信號模型

在不考慮水聲信道時散和多普勒效應的情況下，淺剖回波信號可表示為[11]

其中，s(t)為探測信號，hT(t)為發(fā)射換能器沖激響應，hs(t)為海底淺地層沖激響應，hR(t)為接收換能器沖激響應，n(t)為高斯白噪聲。通常換能器均經(jīng)過校準，因此可將式(1)簡化表示為探測信號與沉積層沖激響應卷積的形式：式(2)中M為沉積層反射界面的個數(shù)am和τm分別為第m反射界面的回波信號幅度和時延。

3 基于FrFT的淺剖精細探測算法

3.1 u域解卷積原理

信號的FrFT可理解為時域信號在時頻平面上逆時針旋轉(zhuǎn)α角度的u域表示，當旋轉(zhuǎn)至αopt= arctan(Δk)+π/2時, Chirp信號會變?yōu)闆_擊信號而出現(xiàn)最佳能量聚集，此時的FrFT稱為最優(yōu)階FrFT，變換階數(shù)popt=2/π稱為最優(yōu)變換階數(shù)[12]。FrFT以量綱歸一化理論為基礎建立算法模型，參量均為量綱歸一化量(無量綱量)，如?k=?B/?T=kS2, ?B, ?T 和S=sqrt(T/fs)分別為量綱歸一化的頻率變化率、帶寬、脈沖寬度和時間因子，符號“?”表示量綱歸一化。

由FrFT公式及其運算性質(zhì)[13,14]，推得Chirp信號關于角度的FrFT為

式中A為常量。由FrFT的階數(shù)可加性，最優(yōu)階FrFT可分解為αc和π/2兩個角度的順序逆時針旋轉(zhuǎn)，而π/2角度坐標軸的逆時針旋轉(zhuǎn)等價于傅里葉變換[15]，因此信號s(t)的最優(yōu)階FrFT又可表示為Sαc(u)的傅里葉變換：

經(jīng)最優(yōu)階FrFT變換，M個淺剖回波信號轉(zhuǎn)變?yōu)镸個u域上的沖擊信號，實現(xiàn)了淺剖信號的u域解卷積處理，得到了沉積層沖激響應的u域表示(下文稱u域沉積層沖激響應)，u域解卷積原理如圖1所示，T為時間窗長度，Ui為Chirp信號的u域分布位置。

實際信號處理時，時間窗長度T是有限的，因此式(5)的積分區(qū)間實際為[-?T/2, ?T/2]，從而淺剖信號最優(yōu)階FrFT應表示為sinc函數(shù)的形式：

圖1 u域解卷積原理

3.2 u域加窗濾波處理

由于Chirp信號為寬帶信號，因此帶內(nèi)信號與噪聲在時頻域上是完全耦合的，采用傳統(tǒng)帶通濾波方法不能消除帶內(nèi)噪聲干擾，而淺剖信號經(jīng)最優(yōu)階FrFT后信號能量主要集中在sinc信號主瓣范圍內(nèi)，高斯白噪聲則分布于幾乎整個u域上，這一過程將耦合區(qū)域由整個頻帶縮小至sinc信號主瓣寬度以內(nèi)，有效實現(xiàn)了信號與噪聲的去耦。從而以u域淺地層沖激響應峰值位置為中心設計窄帶濾波器H(u)，可有效濾除帶內(nèi)噪聲干擾。FrFT為線性變換，因此u域濾波過程可表述為

其中H(u)為幅度歸一化u域濾波器響應。以NH表示濾波器寬度(通常取奇數(shù))，則濾波輸出信噪比為

經(jīng)分析發(fā)現(xiàn)，選擇合適的濾波器窗函數(shù)在有效濾除帶內(nèi)噪聲、最大化輸出信噪比的同時可對沖激響應的旁瓣進行有效抑制，降低高能量信號旁瓣對低能量信號的遮蔽影響，提高弱信號檢測能力。首先，針對4種典型窗函數(shù)進行研究，分析相同窗寬下不同窗函數(shù)對噪聲的抑制效果，見圖2。由于具有寬主瓣和低旁瓣等特性，布萊克曼(blackman)窗噪聲抑制效果最好，較矩形窗輸出信噪比可提高3.8 dB，后續(xù)將采用blackman窗設計u域濾波器。

另外，窗寬是影響沖激響應旁瓣抑制效果的另一個主要因素，不同寬度的blackman窗對旁瓣的抑制效果如圖3所示。對比發(fā)現(xiàn)，當窗寬與sinc信號主瓣寬度相當時(主瓣寬度為17個樣本點)，可抑制掉幾乎所有的旁瓣，同時保留了絕大部分的信號能量。

3.3 淺剖包跡獲取

時頻平面α角度的逆時針旋轉(zhuǎn)可將N點時域Chirp信號變?yōu)镹點u域信號，因此u域信號各樣本點與時域信號各樣本點間必然存在一一映射關系，從而u域樣本序列可通過某一變換關系得到與之對應的時域樣本序列。從量綱歸一化理論出發(fā)，將無量綱u值至時域時間量綱量的變換稱為時間量綱化變換。淺剖包跡是時域沉積層沖激響應的包絡，在獲得高信噪比u域沉積層沖激響應的前提下，通過時間量綱化變換可實現(xiàn)高信噪比淺剖包跡的獲取，從而最終得到高信噪比的淺剖圖像。首先，由式(6)得到淺剖信號幅度可表示為

即u域沉積層沖激響應幅度與時域沉積層沖激響應的幅度間存在一一對應的等比例關系，比例系數(shù)為常數(shù)R。其次，由量綱歸一化理論及Chirp信號在時頻域和u域分布的幾何關系，推得時延為τm的沉積層回波信號在u域的分布位置(或沖激響應峰值位置)為式中fL為Chirp信號起始頻率。為了便于分析，以u軸起始點為零點，則u域各樣本點對應的索引號為

圖3 濾波器寬度對旁瓣的抑制效果

圖2 不同窗函數(shù)時u域濾波輸出SNR

其中，n∈[1,2,…,N]。上式給出了u域樣本序列與時域樣本序列間的映射關系，該式實現(xiàn)了u值的時間量綱化變換。

由圖1幾何關系和式(12)易知，時間量綱化引入了一個時間偏移量ζoffset，該偏移量產(chǎn)生原理如圖4所示，圖中粗虛線表示0時延Chirp信號，其u域分布位置與時域時延相對應，而u軸上大于該位置的區(qū)域與時域 [-ζoffset, 0)相對應。由式(12)也可推得ζoffset為

由式(13)易知，時間偏移量僅與Chirp信號的下限頻率、調(diào)斜率、時間窗長度以及采樣頻率有關，與信號時延無關。因此，時間偏移量可直接計算得到而不受時延信息的影響。

從而時間量綱化得到時域樣本點無偏移的時延值應為

圖4 時間偏移量產(chǎn)生原理

由式(14)可得到u域相鄰樣本點在時域上間隔的樣本點數(shù)為

其中Tp為發(fā)射信號脈寬。當ΔN為整數(shù)時，時間量綱化得到的時域樣本序列具有“插值”的特性(ΔN為插值點數(shù))，因此，通過對該序列做ΔN倍“抽取”和低通濾波即可得到淺剖包跡。綜上，求解淺剖包跡|hs(n)|的過程可表述為

其中，EBPF(n)為帶通濾波輸出回波信號，“Fα”表示p=2α/ π階FrFT“，∠”表示u值的時間量綱化，“||ΔN”表示ΔN倍抽取，LPF表示低通濾波，ceil(N/ΔN)為向趨近于正無窮取整。基于這一求解過程提出了基于FrFT的Chirp淺剖精細探測算法，算法流程如圖5所示。由于淺剖探測為主動聲吶探測，通常已知探測信號參數(shù)，但當考慮多普勒效應時，也可對實際回波信號參數(shù)進行估計以用于算法計算；同時，通過u域沖激響應包絡峰值搜索還可對淺地層回波信號的幅度和時延信息進行估計。

4 算法仿真及數(shù)據(jù)處理分析

4.1 計算機仿真實驗

圖5 淺剖精細探測算法流程圖

采用文獻[16]提出的非等厚沉積層分層模型(層數(shù)可是任意值)，利用式(2)設計淺剖回波信號，通過與脈沖壓縮法和文獻[17]提出的自回歸(AutoRegressive, AR)預測濾波方法進行比對，對該算法性能進行驗證。這里取11層不同厚度的沉積層，由于各層反射強度具有隨機性，回波信號幅度均勻隨機產(chǎn)生，相應的回波時延及歸一化反射強度分別為τ=[3.0, 7.0, 9.5, 11.8, 14.2, 16.0, 18.0, 21.0, 23.0, 25.5, 28.0] ms, a=[1.00, -0.42, 0.67, 0.33, -0.50, 0.25, 0.34, 0.28, -0.50, -0.37, 0.25]。發(fā)射信號為脈寬6 ms、頻率2～6 kHz的Chirp信號，信噪比6 dB，觀測時間窗長度36 ms。由式(15)計算時間量綱化后采用2倍抽取(ΔN=2)，則fs=41.56 kHz。

實驗1 將帶通濾波輸出淺剖信號做最優(yōu)階(1.33318階)FrFT變換，得到u域沉積層沖激響應并進行u域加窗濾波處理，blackman窗寬為15點，濾波效果如圖6所示，濾波結(jié)果表明：u域解卷積和加窗濾波在有效濾除帶內(nèi)噪聲的同時對沖激響應旁瓣進行了有效抑制。由式(13)計算時間偏移量為15.025 ms，對圖6中高信噪比的u域沉積層沖激響應包絡進行時間量綱化變換，經(jīng)2倍抽取和低通濾波得到了高信噪比的淺剖包跡信號，如圖7所示，與脈沖壓縮方法相比，本文算法輸出時域波形具有更低的旁瓣和更高的信噪比，且精確探測到了不同反射強度的沉積層。

同時，采用AR預測濾波方法對同一組信號進行處理，處理結(jié)果見圖8，由于AR預測濾波基于傳統(tǒng)的解卷積方法，對噪聲干擾較敏感，雖能探測到主要的淺地層信息，但帶內(nèi)噪聲干擾仍較嚴重。仿真實驗從u域和時域上均表明本文算法對帶內(nèi)噪聲和沖激響應旁瓣的抑制效果，直觀地反映了其優(yōu)于脈沖壓縮和AR預測濾波方法的淺剖精細探測能力。

圖7 淺剖包跡信號及其頻譜

圖8 AR預測濾波輸出淺剖包跡及其頻譜

實驗2 輸入信噪比為-10～10 dB時，采用Monte Carlo方法算法輸出信噪比進行5000次統(tǒng)計實驗，實驗結(jié)果見圖9。由于有效抑制了帶內(nèi)噪聲干擾，本文算法輸出信噪比最高，與脈沖壓縮法相比，在低于0 dB的輸入信噪比時可獲得10 dB左右的增益。帶內(nèi)噪聲的有效抑制有利于對弱回波信號的檢測，從而有利于獲得更精細的淺地層剖面圖像。

圖9 輸出信噪比變化曲線

4.2 數(shù)據(jù)處理分析

采用兩組Chirp淺剖原始數(shù)據(jù)對算法性能進行驗證，Chirp探測信號頻率范圍2～6 kHz，脈寬32 ms，采樣頻率26.316 kHz，處理結(jié)果如圖10和圖11所示。

圖10是存在斜坡和溝壑等復雜地形的沉積層剖面處理結(jié)果。結(jié)果表明，本文算法可以有效地檢測到斜坡地形的弱回波信號(圖中實線框區(qū)域)，雖然AR預測濾波法也檢測到了弱回波信號，但背景噪聲明顯高于本文算法，淺剖圖像圖10(c)較圖10(a)和圖10(b)更加清晰，細微地形更加豐富，充分說明了算法對帶內(nèi)噪聲較強的抑制能力。

圖11是對多層沉積層的探測結(jié)果，脈沖壓縮法和AR預測濾波法僅能探測到海底以下20 m左右的地層結(jié)構，無法探測更深的沉積層結(jié)構特征。圖11(c)中不僅清晰地探測到了海底以下20 m區(qū)域的沉積層，而且還探測到了100 m深度的弱回波信號，同時整個探測結(jié)果的背景噪聲干擾抑制效果明顯，淺剖圖像更加清晰，計算其信噪比比圖11(a)和圖11(b)分別提高了約7.2 dB和8.4 dB。為了更加直觀地說明本文算法處理性能，隨機選取第495個探測周期的數(shù)據(jù)采用上述3種方法進行處理，處理結(jié)果如圖12至圖14所示，時域波形及淺剖包跡處理結(jié)果均表明了本文算法的噪聲抑制效果以及地層的精細探測能力，特別實現(xiàn)了100 m深度處弱回波信號的檢測。數(shù)據(jù)2的處理結(jié)果再次驗證了本文算法優(yōu)于其它兩種方法的處理性能。

圖10 數(shù)據(jù)1處理結(jié)果

圖11 數(shù)據(jù)2處理結(jié)果

圖13 AR預測濾波處理結(jié)果

圖12 脈沖壓縮處理結(jié)果

圖14 本文算法處理結(jié)果

5 結(jié)束語

本文在對FrFT解卷積原理、u域加窗濾波以及時間量綱化理論研究的基礎上，提出了一種基于FrFT的Chirp淺剖精細探測新方法，該方法在u域?qū)崿F(xiàn)海底淺地層的解卷積和濾波處理，具有優(yōu)于傳統(tǒng)方法的帶內(nèi)信號與噪聲的去耦與濾波效果，從而具有更高的信噪比增益以及對微弱地層回波信號的檢測能力，時間量綱化變換將高信噪比的u域沉積層沖激響應包絡變換至時域，直接得到了高信噪比的淺剖包跡，并且由于u域加窗濾波處理，本文算法輸出淺剖包跡較脈沖壓縮輸出沖擊信號具有更強的旁瓣抑制能力。計算機仿真和實驗數(shù)據(jù)處理驗證了該方法優(yōu)于脈沖壓縮和AR預測濾波方法，為淺地層剖面精細探測提供了新的技術途徑。

[1] Gron O and Boldreel L O. Sub-bottom profiling for large-scale maritime archaeological survey[C]. Proceedings of Oceans 2013 Marine Technology Society, Bergen, Norway, 2013: 1-8.

[2] Sandrine R, Michel L, Rozenn D, et al.. Forward modeling for marine sediment characterization using chirp sonars[J]. Geophysics, 2011, 76(4): 91-99.

[3] Tseng Yao-ting, Ding Jian-jun, and Liu Char-shine. Analysis of attenuation measurements in ocean sediments using normal incidence chirp sonar[J]. IEEE Journal of Oceanic Engineering, 2012, 37(3): 533-543.

[4] Xie Yu-jiang and Liu Gao. Predicting the distribution of reservoirs by combining variable wavelet model of seismograms with wavelet edge analysis and modeling[J]. Journal of Applied Acoustics, 2014, 101: 116-123.

[5] Hakan K and Derman D. A procedure to reduce side lobes of reflection wavelets: a contribution to low frequency information[J]. Journal of Applied Geophysics, 2013, 96: 107-118.

[6] 齊林, 陶然, 周思永, 等. 基于分數(shù)階傅里葉變換的線性調(diào)頻信號的自適應時頻濾波[J]. 兵工學報, 2003, 24(4): 499-503. Qi Lin, Tao Ran, Zhou Si-yong, et al.. An adaptive time-frequency filtering method based on fractional Fourier transform for linear frequency modulation signals[J]. Acta Armamentarii, 2003, 24(4): 499-503.

[7] Hakan K. Further improvement of temporal resolution of seismic data by autoregressive (AR) spectral extrapolation[J]. Journal of Applied Geophysics, 2006, 59(4): 324-336.

[8] Zeng Wen-jun, Jiang Xue, Li Xi-lin, et al. Deconvolution of sparse underwater acoustic multipath channel with a large time-delay spread[J]. The Journal of the Acoustical Society of America, 2010, 127(2): 909-919.

[9] 朱健東, 趙擁軍, 唐江. 線性調(diào)頻連續(xù)波信號的周期分數(shù)階Fourier變換檢測與估計[J]. 電子與信息學報, 2013, 35(8): 1827-1833. Zhu Jian-dong, Zhao Yong-jun, and Tang Jiang. Periodic FRFT based detection and estimation for LFMCW signal[J]. Journal of Electronics & Information Technology, 2013, 35(8): 1827-1833.

[10] 陳鵬, 侯朝煥, 馬曉川, 等. 基于匹配濾波和離散分數(shù)階傅里葉變換的水下動目標LFM回波聯(lián)合檢測[J]. 電子與信息學報, 2007, 29(10): 2305-2308. Chen Peng, Hou Chao-huan, Ma Xiao-chuan, et al.. The joint detection to underwater moving target’s LFM echo based on matched filter and discrete fractional fourier transform[J]. Journal of Electronics & Information Technology, 2007, 29(10): 2305-2308.

[11] Carcreff E, Bourguignon S, Idier J, et al.. Including frequency-dependent attenuation for the deconvolution of ultrasonic signals[C]. Proceedings of Meetings on Acoustics, Montreal, Canada, 2013, 19: 1-9.

[12] 仇兆煬, 陳蓉, 汪一鳴. 基于FRFT的線性調(diào)頻信號欠采樣快速檢測方法[J]. 電子學報, 2012, 40(11): 2165-2170. Qiu Zhao-yang, Chen Rong, and Wang Yi-ming. Fast detection of LFM signal based on FRFT and sub-Nyquist sampling[J]. Acta Electronical Sinica, 2012, 40(11): 2165-2170.

[13] 李麗, 邱天爽. 基于分數(shù)階傅里葉變換的雙基地雷達線調(diào)頻信號的參數(shù)聯(lián)合估計新方法[J]. 電子與信息學報, 2012, 34(4): 878-884. Li Li and Qiu Tian-shuang. A novel method for joint parameter estimation of LFM signals in bistatic MIMO radar system on FRFT[J]. Journal of Electronics & Information Technology, 2012, 34(4): 878-884.

[14] Khalid H M and David K P. Geophysical signal parameterization and filtering using the fractional Fourier transform[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(3): 845-852.

[15] 阮航, 吳彥鴻, 葉偉. 基于FRFT-CLEAN的機動目標逆合成孔徑激光雷達成像算法[J]. 電子與信息學報, 2013, 35(7): 1540-1546. Ruan Hang, Wu Yan-hong, and Ye Wei. Inverse synthetic aperture radar imaging algorithm for maneuvering target based on FRFT-CLEAN[J]. Journal of Electronics & Information Technology, 2013, 35(7): 1540-1546.

[16] Quinquis A and Radoi E. FM pulses separation for improving sub-bottom attenuation estimation[C]. Proceedings of Oceans’97, Halifax, Canada, 1997, 2: 1359-1365.

[17] Satchidananda P. Remote acoustic evaluation of seafloor sediment properties[D]. [Ph.D. dissertation], University of Rhode Island, 1992.

朱建軍： 男，1981年生，助教，博士生，研究方向為海洋工程信息系統(tǒng)、水聲信號處理.

魏玉闊： 男，1982年生，講師，博士，研究方向為海洋工程信息系統(tǒng)、水聲工程.

杜偉東： 男，1984年生，博士生，研究方向為海洋工程信息系統(tǒng).

Chirp Sub-bottom Profiling Detailed Detection Method Based on Fractional Fourier Transform

Zhu Jian-jun①②Wei Yu-kuo①②Du Wei-dong②Li Hai-sen①②①( Acoustic Science and Technology Laboratory, Harbin Engineering University, Harbin 150001, China)
②(College of Underwater Acoustic Engineering, Harbin Engineering University, Harbin 150001, China)

Weak signal detection and high SNR seismic image generation are primary tasks in detailed sub-bottom profile detection. After analyzing the principle of deconvolution based on Fractional Fourier Transform (FrFT) and deriving the formula of time dimensional transformation, a new detailed sub-bottom profile detection algorithm based on FrFT is proposed. The fractional Fourier domain (u domain) sub-bottom impulse response is achieved by u domain deconvolution and the intraband SNR is increased by u domain windowed filtering technique, then high SNR envelop of u domain sediment impulse response envelop is transformed to time domain by time dimensional transformation to get high quality sub-bottom profile. Simulation and experimental data processing validate the validity of the algorithm in intraband denoising and detailed detection, and its performance is better than pulse compression and AutoRegressive (AR) forecast filtering.

Signal processing; Sub-bottom profile; Detailed detection; Fractional Fourier Transform (FrFT); Intraband denoising; Time dimensional transform

TN911.7

A

1009-5896(2015)01-0103-07

10.11999/JEIT140140

2014-01-21收到，2014-04-15改回

國家國際科技合作項目(2008DFR70320)，國家自然科學基金(41327004, 41306182, 61401112)， 教 育 部 博 士 點 基 金(20112304130003)和哈爾濱市科技攻關項目(2008AA2AE005)資助課題

*通信作者：李海森 hsenli@126.com