梅曉東,紀永剛,王祎鳴,唐慶樂
(1.中國石油大學(華東),山東 青島266580;2.國家海洋局 第一海洋研究所,山東 青島266061;3.國家海洋局 航天科技集團公司海洋遙測工程技術研究中心,山東 青島266061)
海洋表面流的流速是根據回波多普勒譜中一階峰相對于布拉格頻率的偏移量計算的[1]。窄波束地波雷達一階峰是很窄的尖峰,因此分離其一階譜區可通過直接搜索一階峰位置或利用參數譜估計方法定位實現[2-3]。而寬波束雷達向海面發射電磁波,接收同一時刻來自不同方位的回波信號,當不同方位海流流速不同時,回波譜中各一階峰會出現不同偏移,導致一階譜區呈現出展寬特性[4]。因此利用寬波束地波雷達提取流速,分離一階譜區是重要的步驟,它影響著各距離元的流速范圍及流速分布結果的準確性。現有的基于質心估計法和差譜法的一階譜區分離方法存在穩定性差的問題,質心估計法[5-6]簡單直接,由于一階譜寬度設為固定經驗值,自適應性較差;差譜法[7]在一階峰之間幅度變化較大時效果不理想。
為了獲取遠距離海態信息,雷達通常選取較低的發射頻率,因此接收到的散射回波譜中二階分量易被背景噪聲掩蓋,一階譜幅度顯著高于周圍較為平坦的噪聲。基于這一特點,我們提出一種信噪比方法來提取寬波束地波雷達的一階回波譜區,為地波雷達海流反演提供準確可靠的一階回波譜信息。基于實測地波雷達數據的處理結果及與差譜法的對比分析,得出該方法的一般性能,最后針對信噪比參數的選擇問題進行了簡要分析研究。
低頻段的高頻地波雷達回波多普勒譜示意圖見圖1。通過觀察可發現2個一階分量分別位于零階峰兩側且幅度遠高于平均水平,回波譜中并未出現明顯的二階譜,總體來看一階峰以外的連續譜較為平坦,波動較小。當在圖1所示多普勒譜上某位置截取一定寬度多普勒譜作為信號,相鄰一定寬度的連續譜作為噪聲,利用這2個參數求取信噪比,在遍歷多普勒譜所有值后,信噪比的極大值一定是出現在一階譜位置上。基于信噪比的一階譜區分離方法具體流程詳見圖2。我們利用高頻地波雷達左右兩個一階回波譜中心間距不變的特性,根據實際最大海流流速引起一階回波譜頻移及展寬大小來確定信號窗及噪聲框的取值范圍,并依此確定一階譜區信號窗中心位置和信號窗口的取值范圍;然后滑動信號窗和噪聲窗位置計算一階回波譜信號和噪聲的信噪比,在此過程找出的信噪比的極大值,即可確定此時信號窗和噪聲窗的位置,最終確定一階譜的范圍。

圖1 信號窗與噪聲窗的示意圖Fig.1 A sketch of signal windows and noise windows

圖2 信噪比法提取一階譜流程圖Fig.2 Procedures of the separation of first-order spectra with the SNR method
設回波多普勒譜為Doppler(n),n=1,2,…,N。為求取信噪比,首先將信號部分用寬度為l的矩形窗進行框定(圖2中的實線矩形窗),窗體中心位置設為a,信號窗內的多普勒幅度為所需信號。滑動窗長度為l。設多普勒分辨率為δf,布拉格頻率為fB,最大可能流速引起的頻率偏移量為Δf,那么左側一階譜區信號窗中心位置a的取值范圍設為a1~a2,其中,a1=-fB-Δf/2,a2=-fB+Δf/2。信號窗寬度l的頻率取值范圍設為l1~l2,其中l1=2δf,l2=Δf。設噪聲窗與信號窗的寬度比值為k,則噪聲窗寬度為k·l,圖1中,黑色矩形窗表示信號窗,其寬度是l,表示占據了l個多普勒點,虛線所標注的是噪聲窗,其寬度是l的k倍。
由于多普勒譜在左右兩側有2個對稱的一階譜,因此需要同時計算。首先,左側一階譜區信號可表示為

式中b1=a-0.5l,b2=a+0.5l。如圖2所示,與信號窗緊鄰的2個虛線矩形表示噪聲的范圍,寬度均為k·l。左一階譜區周圍噪聲的表達式可寫為

式中,c1=a-(k+0.5)l,c2=a+(k+0.5)l。利用左右一階峰間距為2fB的特性,右側一階譜區信號窗中心位置取值范圍為(fB-Δf/2)~(fB+Δf/2)。其信號和噪聲可表示為

其中,d1=a+2fB-0.5l,d2=a+2fB+0.5l;e1=a+2fB-(k+0.5)l,e2=a+2fB+(k+0.5)l;f2=a+2fB-0.5l,f2=a+2fB+0.5l。當左右一階譜區的信號窗與噪聲窗在多普勒譜上同時進行滑動,即參數a與l分別在a1~a2和l1~l2范圍內變化時,利用式(2)和(3)可分別得到2組信噪比值SNRR和SNRL。為了保證兩側一階譜區結果的一致性,定義最終信噪比為

選擇其中的最大值max(SNR)?(aSNRmax,lSNRmax),對應的aSNRmax和aSNRmax+2fB分別為左、右一階譜區的中心位置,窗長l為一階譜寬度,利用周圍噪聲與一階譜的對應關系(k·l)可求出噪聲范圍,進而計算出確定譜峰所需閾值。只有位于一階譜邊界內,且幅度大于閾值的多普勒信號可用于流速的提取。
實測數據來自威海的某陣列式高頻地波雷達數據,雷達主要探測海域在渤海,雷達頻率為8.9MHz,8陣元,陣元間距14.5m,數據積累時間約為9min。多普勒分辨率為0.002Hz,布拉格頻率為0.304Hz,最大可能流速為150cm/s,對應的頻偏為0.089Hz。左側一階譜區信號窗中心位置a頻率范圍為-0.349~-0.260Hz,信號窗寬度l的頻率取值范圍為0.004~0.089Hz。取k=0.5,噪聲窗寬度為0.002~0.044 Hz。右側一階譜區參數可利用左右一階峰頻率間隔為2fB的特性推得。
為驗證信噪比方法的有效性,利用該方法對實測的高頻地波雷達數據進行處理,并將該方法的處理結果與差譜法的處理結果進行對比分析。差譜法原理[7]為將各通道合成的回波多普勒功率譜采用3點滑動平均,可在保持回波波譜分布的前提下,減少可能出現的一階峰過度分裂以及緊鄰的高階譜毛刺造成的后續過程的判斷失誤;取滑動平均后的功率譜的對數,得到對數譜;將對數譜相鄰點作差,得到差譜;將差譜作3點滑動平均;利用一階譜區相對于緊鄰的高階譜區有較陡峰的特征,取差譜中正負布拉格頻率附近的最大正值對應于一階譜區的左邊界,最小負值對應于一階譜區的右邊界。分別計算分離出的左右一階譜所夾區域以外的平均電平作為左右區域的局部噪聲電平,將左右局部噪聲電平加上一定信噪比,作為一階譜區內判定可用信號的局部閾值,局部閾值以上的為可提取流速的信號。3種不同距離元數據的處理結果見圖3、圖4和圖5。實線矩形框表示信噪比方法的一階譜區分離結果,虛線框表示差譜法一階譜區分離結果。
分析發現:1)大部分情況下,信噪比方法和差譜法確定的一階譜邊界范圍差別很小,分離的一階峰相同。第45個距離元的回波信號處理結果見圖3,左右一階譜區目視范圍分別為-0.347~-0.291和0.258~0.314Hz(圖3c)。由局部放大圖3c可以看出,信噪比方法(實線)所確定左側一階譜范圍為-0.341~-0.310Hz,而差譜法是-0.332~-0.316Hz,邊界雖不完全相等,但其中所包含一階峰數目都為3個,對應徑向流速分別為-43.99cm/s,-31.52cm/s,-25.28cm/s,以遠離雷達站方向運動。因此最終的流速個數并無差異;右側的處理結果與左側類似。由此判定,2種方法所確定的邊界范圍與目視范圍吻合度較好,未出現流速漏檢或誤判。
2)海洋表面具有隨機波動性,因此不同方位上的海流流速可能表現出較大差異,在多普勒譜上表現為一階峰之間的明顯分離。左、右一階譜區為-0.348~-0.293Hz和0.264~0.304Hz,但在左側一階譜區的-0.328Hz位置處一階峰值之間出現分裂(圖4)。信噪比方法(實線)所確定的一階譜區范圍是-0.338~-0.314Hz,而差譜法受一階峰影響,右邊界恰好位于分裂處(圖4c)。這是由于,差譜法所求一階譜區右邊界應位于差譜中布拉格頻率附近最小值。分析差譜圖(圖6)可發現,最小值剛好位于-0.330Hz點,即一階峰分裂點。這就導致差譜法出現判別錯誤,邊界未能準確覆蓋所有一階峰,發生流速漏檢情況。圖4d給出的是右一階譜區的放大圖,可以看到2種方法得到的邊界是重合的,這是由于右側的一階峰之間較為連續,無大的分裂,表明2種方法的判斷結果較為準確和一致。

圖3 第45距離元回波多普勒及一階譜區分離結果Fig.3 The Doppler spectra of echoes from No.45range cell and the separation results of the first-order spectral region

圖4 多普勒譜內部出現分裂時的處理結果Fig.4 Processing results when split occurs within the Doppler spectra

圖5 一階譜附近出現硬目標回波情況下的處理結果Fig.5 Processing results when hard target echo occurs around the first-order spectra

圖6 差譜圖Fig.6 Differential spectra
3)當海上存在艦船目標時,雷達接收的后向散射回波不僅包含海面狀態信息,多普勒譜中還會出現硬目標回波,對流速的提取產生一定干擾。在緊鄰右側一階譜區的地方出現一幅度較強的回波信號,頻率位置為0.317Hz(圖5a)。差譜法受到該回波信號的影響,將其誤判為流速。圖7中硬目標回波處的差譜值明顯小于目視邊界點,可以解釋差譜法出現誤判的原因。而信噪比法確定的一階譜區邊界與目視范圍較為接近,且未受到硬目標回波的干擾,效果較好。

圖7 差譜圖Fig.7 Differential spectra
我們提出的信噪比方法利用一階譜信號與周圍噪聲之間信噪比最大的特性,式(1)~(6)對該方法的原理及實現過程進行描述,其中噪聲寬度為信號寬度的k倍。為分析噪聲窗寬度不同對處理結果的影響,這里分別在k=0.2,k=0.5和k=1.2時對同一組實測高頻地波雷達數據進行了處理(圖8)。從處理結果來看,k取不同值對應的一階譜邊界位置雖不完全相同,但差異不大,甚至出現重合,且由于各邊界劃定范圍內的一階峰個數相同,因此未對流速結果造成影響。

圖8 k值對階譜區邊界確定的影響Fig.8 The influence of kvalue on the determination of the boundary of the first-order spectral region
基于低頻段地波雷達的海洋回波譜特性,提出一種新的一階譜區邊界確定方法;選取3組典型實測雷達回波數據,利用信噪比法和差譜法分別進行處理,從對比2種處理結果可以看出,我們提出的方法可準確分離一階譜區,并避免硬目標回波和一階峰之間分裂較大的影響;通過對比噪聲窗寬度取不同值對應的一階譜區邊界,發現噪聲窗與信號窗的寬度比值k對分離結果的影響很小。
在應用于發射頻率較低、二階譜不明顯的地波雷達回波多普勒譜時,與差譜法相比,我們的方法優勢突出,具有較好的處理效果。但當雷達頻率較高時,會激起較高幅度的二階譜,此時一階譜與周圍連續譜之間信噪比降低,在利用新方法時可能出現誤判情況,因此需要進一步發展相應處理方法。
(References):
[1]ZHOU Z X,LIU Y T.Computational method for the extraction of surface-current by HF ground-wave-radar[J].Marine Science Bulletin,1997,16(5):79-84.周志鑫,劉永坦.高頻地波雷達提取表層海流的計算方法[J].海洋通報,1997,16(5):79-84.
[2]HICKEY K,KHAN R,WALSH J.Parametric estimation of ocean surface currents with HF radar[J].IEEE Journal of Oceanic Engineering,1995,20(2):139-144.
[3]MARTIN R J,KEARNEY M J.Remote sea current sensing using HF radar:an autoregressive approach[J].IEEE Journal of Oceanic Engineering,1997,22(1):151-155.
[4]ZHU D Y.Applications of high frequency ground wave radar to coastal ocean[D].Xiamen:Xiamen University.2008.朱大勇.高頻地波雷達在近海區域的應用研究[D].廈門:廈門大學,2008.
[5]LEISE J A.The analysis and digital signal processing of NOAA's surface current mapping system[J].IEEE Journal of Oceanic Engineering,1984,9(2):106-113.
[6]MILLER P A,LEISE J A.Radar Doppler detection methods with applications to CODAR[M].[S.l.]:NOAA,Environmental Research Laboratories,1982:25.
[7]YANG S L,KE H Y,HOU J C,et al.Signal preprocessing for bearing determination of ocean surface radial current mapping based on MUSIC[J].Modern Radar,2001,23(4):51-52.楊紹麟,柯亨玉,侯杰昌,等.MUSIC算法提取海洋表面徑向流方位的信號預處理[J].現代雷達,2001,23(4):51-52.