利用靜態(tài)方向圖數(shù)值綜合方法抑制電離層干擾

陳驍鋒 吳雄斌 徐興安 沈志奔 王 立

(武漢大學電子信息學院海態(tài)實驗室，湖北 武漢 430079)

引 言

電離層干擾一直是高頻地波雷達探測海洋過程中非常重要的研究課題. 國內(nèi)外關于此方面的研究主要圍繞兩個方面：一是電離層雜波特性的研究[1]；二是電離層雜波的抑制方法研究[2]. 雖然國內(nèi)外已經(jīng)提出了許多抑制電離層干擾的方法，但由于電離層干擾的復雜性，大多數(shù)文獻都是針對某段特殊的干擾數(shù)據(jù)進行分析，具有很強的針對性. 文獻[3-5]嘗試從陣列信號處理入手，通過控制雷達波束的方式，限制來自空中的電離層干擾，達到抑制的目的[3-5].

自適應陣列處理和空間譜估計是陣列信號處理的兩個主要研究方向. 空間譜估計算法能夠突破“瑞利限”實現(xiàn)高分辨的到達角方位估計，高頻地波雷達在探測海洋表面流時，采用的就是空間譜估計算法中的多重信號分類法(Multiple Signal Characteristic,MUSIC). 然而在某些特定的環(huán)境下，空間譜估計算法由于存在探測信號源數(shù)必須小于天線陣元數(shù)的要求[4]，會導致算法的失效，如利用寬波束雷達探測海洋表面的海浪. 此類情況下必須考慮用波束形成的方法，形成指定方向較窄波束的回波譜[5].

在理想條件下，自適應波束形成技術可以有效抑制干擾，使陣列輸出信干噪比最大. 但實際系統(tǒng)中，不可避免的存在各種誤差，如陣列通道的幅相誤差、互耦誤差、位置誤差等，自適應的波束形成對這些誤差干擾非常敏感[6]. 靜態(tài)波束形成的方法與之相比，在實際應用中更加穩(wěn)定，具有較好的魯棒性.

空間匹配濾波器又稱常規(guī)波束形成(Conventional Beamforming，CBF)是最基礎的靜態(tài)波束形成方法，它用對各陣元信號加權進行相位補償?shù)姆椒ㄊ怪鞑ㄊ赶蚱谕男盘柗较? Dolph于1946年提出了用于均勻線陣的經(jīng)典低副瓣綜合方法(Dolph-Chebyshev方法)，該方法可以在限定最大副瓣電平的同時獲取較窄的主瓣波束[7]. 1966年Tseng和Baklanov將Dolph-Chebyshev方法推廣應用于均勻矩形陣列[8-9]. 1990年Olen提出一種較為簡單、直觀的基于自適應陣列理論的靜態(tài)方向圖數(shù)值綜合方法(Numerical Pattern Synthesis ，NPS[10]). 1992年Er提出的峰值控制方法可以在任意陣列情況下進行復雜靜態(tài)方向圖綜合[11]，但由于在每次迭代中需確定方向圖的峰值點，因此較Olen的方法復雜. 2008年Eom利用遺傳算法計算信號的最大功率，對相控陣天線單元的相位和幅度進行了優(yōu)化和校正[12]，2010年Frank用遺傳算法形成二維平面陣的窄波束平覆蓋區(qū)域[13].

在實際應用過程中，高頻地波雷達回波信號處理中的波束形成方法，都只考慮一維即水平方向上的影響，這樣導致的結(jié)果是盡管水平方向圖具有窄波束、低旁瓣等優(yōu)勢，但來自空中的干擾也會大量進入到回波信號中，其中最典型的就是電離層干擾[14]. 特別對于小陣列而言，因為場地有限，往往只能采用雙排陣列，導致陣列方向圖在俯仰角方向上的分辨率非常低，嚴重影響回波特性.

針對這一問題，本文提出利用二維NPS算法對高頻地波雷達回波信號進行波束形成，并針對NPS迭代過程中的性能變化特性，提出一種基于方向圖旁瓣誤差的性能模型，能夠用于尋找NPS的最優(yōu)設置參數(shù). 分別針對方陣和特殊小型陣列進行了仿真分析，通過實測數(shù)據(jù)的比較，能夠發(fā)現(xiàn)改進的二維NPS算法能夠有效抑制來自垂直方向的干擾，為海洋信息的提取提供可靠的保障.

1 原理算法

首先介紹任意陣列雷達信號的數(shù)學建模，分析利用二維NPS算法進行波束形成過程中的主要步驟及關鍵問題.

對于一個N元任意窄帶信號，導向矢量可以表示為

a(θ,φ)=[g1(θ)ejφ1(θ,φ),g2(θ)ejφ2(θ,φ),

…,gN(θ)ejφN(θ,φ)].

(1)

式中：gi(θ)為第i個陣元的響應；θ表示水平入射方位角；φ表示俯仰角；φi(θ,φ)為第i個陣元相對相位參考點的相移. 設陣列權矢量W=[w1,w2, …,wN]T， 則陣列響應方向圖可表示為

G(θ,φ)=|WHa(θ,φ)|.

(2)

設(θ0,φ0)為期望信號角度，即歸一化處理后有G(θ0,φ0)=1. 第一次迭代時的初值為

(3)

設D(θ,φ)為期望的旁瓣峰值電平，以dB表示，則實際電壓值為

d(θ,φ)=10D(θ,φ)/20.

(4)

根據(jù)經(jīng)典的NPS算法，在陣列的旁瓣區(qū)施加多個虛擬干擾，并控制虛擬干擾功率的大小，經(jīng)過迭代運算，適當?shù)恼{(diào)整該位置的干擾強度，從而得到期望的陣列天線方向圖. 若第k次迭代求得的自適應權矢量為WH(k)，則第k次迭代的方向圖增益為

G(θ,φ,k)=|WH(k)a(θ,φ)|

(5)

根據(jù)式(5)，可以確定出主瓣范圍. 由于二維方向圖的主瓣是一個放射性的波束，主瓣范圍的確定相對一維來說復雜許多. 雖然一些方法如圖像處理技術中的灰度閾值法、模式識別法[15]等，能夠更準確地確定主瓣范圍，但此時主瓣范圍的表示形式也更加復雜. 因為每次迭代都需要進行主瓣區(qū)域的確定，所以，需要根據(jù)實際情況采用合適的方法，使時間復雜度不至于過大，又保證迭代結(jié)果的質(zhì)量.

為了算法簡便，只提取陣列水平平面和法線方向的垂面這兩個平面，利用差分法分別檢測主瓣附近的變化趨勢，從而確定包絡主瓣的方形區(qū)域[θL,θR,φL,φR].

假設M個干擾均勻分布在整個空間區(qū)域，入射角掃描個數(shù)為P，俯仰角掃描個數(shù)為Q，M=P×Q. 對于處于主瓣范圍內(nèi)的虛擬干擾，令其功率為0；對于主瓣范圍外的干擾，通過比較虛擬干擾環(huán)境對應的方向圖在旁瓣(θm,φm)處的增益G(θm,φm)與期望的旁瓣增益d(θm,φm)的大小，對虛擬干擾的功率進行調(diào)整. 虛擬干擾功率的迭代公式為

Γ(θp,φq,k)=ξ(θp,φq,k)+η[G(θp,φq,k)-

d((θp,φq))]，

ξ(θp,φq,k+1)=

(6)

式中:η為一合適的常數(shù)因子，η的取值既要保證迭代的收斂性，又要保證收斂的速度，通常通過試驗經(jīng)驗獲得;p和q表示方位.

因為假設的虛擬干擾和噪聲互不相關，可以得到第k次迭代的陣列協(xié)方差矩陣

aH(θp,φq)].

(7)

(8)

最終迭代結(jié)束時得到的自適應權矢量即為所求的靜態(tài)權矢量.

2 仿真與實測數(shù)據(jù)分析

2.1 性能仿真

為驗證二維NPS波束形成算法的性能及其對OSMAR071天線陣列的適用性，利用兩種不同大小、規(guī)格的陣型(8×8方陣和OSMAR071天線陣)對算法進行仿真分析，給出兩種陣列的二維CBF和NPS的方向圖結(jié)果以便比較. 同時仿真不同方向圖下存在電離層干擾的信號反演結(jié)果，利用蒙特卡洛算法比較兩種方法的性能.

1) 8×8元方陣

取8×8元方陣，陣元間距為半波長，主瓣指向陣列法線方向(方位角90°，俯仰角0°)，設置旁瓣峰值為-40 dB，迭代次數(shù)k=1 000次，每隔一個方位角、一個俯仰角取一個虛擬干擾，迭代的常數(shù)因子為η=0.1.

圖1(a)、(b)所示為8×8元陣列二維CBF和二維NPS仿真結(jié)果的比較. 8×8元陣列二維NPS的結(jié)果與二維CBF的結(jié)果相比：①NPS水平方向的旁瓣波束得到有效的降低，旁瓣電平下降了約7 dB；②NPS俯仰角方向上的旁瓣得到有效抑制，特別對于俯仰角較大的區(qū)域，能夠達到10 dB以上的抑制.

(a) 8×8元陣二維CBF方向圖

(b) 8×8元陣二維NPS方向圖

(c) OSMAR071二維CBF方向圖

(d) OSMAR071二維NPS方向圖圖1 兩陣列二維CBF和二維NPS仿真結(jié)果

2) OSMAR071天線陣

OSMAR071天線陣的陣型坐標如表1所示，工作頻率為8 MHz，主瓣指向陣列法線方向(方位角90°，俯仰角0°)，設置旁瓣峰值為：-15 dB(方位角0°～180°，俯仰角0°～5°)，-10 dB(方位角0°～180°，俯仰角6°～90°)，迭代次數(shù)k=1 000次，每隔一個方位角、一個俯仰角取一個虛擬干擾，迭代的常數(shù)因子為η=0.1.

表1 OSMAR071天線陣天線位置示意圖 m

圖1(c)、(d)所示為OSMAR071天線陣二維CBF和NPS仿真結(jié)果的比較. OSMAR071天線陣二維NPS的結(jié)果與二維CBF結(jié)果相比：① NPS水平方向的波束旁瓣電平相比CBF要高2～3 dB，且主瓣寬度明顯增大，達到26°，-3 dB處主瓣寬度變大2°；②NPS俯仰角方向上的旁瓣得到有效抑制，能夠達到5 dB以上. 產(chǎn)生這一現(xiàn)象的原因是OSMAR071的小型天線陣列只有兩排，在空間的上分辨率非常低，水平方向和垂直方向的方向圖相互制約，強制優(yōu)化俯仰角方向的旁瓣，會導致水平方向的方向圖優(yōu)化效果變差.

由圖1可以發(fā)現(xiàn)：兩種算法在不同陣列下的效果是不同的，在8×8元陣列中，二維NPS算法能夠很容易實現(xiàn)方向圖的優(yōu)化；而在OSMAR071小型陣列中，在對俯仰角較大區(qū)域的干擾進行抑制的時候，由于小型陣列的陣型限制，不可避免地會使水平方向的性能下降.

在仿真的過程中還發(fā)現(xiàn)，通過控制不同方位的虛擬干擾功率大小，得到的方向圖效果是不同的，仿真設置是虛擬功率為-15 dB(方位角0°～180°，俯仰角0°～5°)和-10 dB(方位角0°～180°，俯仰角6°～90°)，只是一個粗糙的仿真設置. 如果合理選用一些參數(shù)優(yōu)化算法，對虛擬干擾平面進行細致的設置，會得到性能更好的方向圖.

同時在迭代過程中，方向圖的變化并不是平緩的，在8×8元陣列k=504及k=918仿真時，會出現(xiàn)旁瓣電平突然變大，然后再緩慢減小的現(xiàn)象.

這些現(xiàn)象說明，在利用二維NPS對陣列進行波束形成時，必須針對實際情況，設置更加合理的參數(shù)設置，以尋求最合適的優(yōu)化結(jié)果.

3) 電離層干擾抑制性能仿真

為了進一步比較兩種算法對電離層干擾抑制的效果，特構造一個存在噪聲時的正弦信號模型. 該信號的幅度為1，頻率為歸一化-1，方位角為60°，俯仰角為0° ；電離層干擾的幅度為0.000 1，頻率為覆蓋頻段上的均勻隨機分布，方位角覆蓋寬度為60°，俯仰角為80°；隨機白噪聲的幅度為0.001.

圖2(b)中的電離層干擾的方位角位置為-30°～30°，圖2(c)中的電離層干擾方位角為-180°～180°內(nèi)隨機起始位置，覆蓋寬度為60°. 可以看出二維NPS算法與二維CBF算法相比，能夠更好的抑制來自空中各個方向的電離層干擾，提高信號的信噪比.

(a) 無電離層干擾下波束形成圖

(b) 電離層干擾下波束形成圖

(c) 1 000次實驗噪聲變化圖2 電離層干擾抑制性能仿真

針對參數(shù)選擇的問題，提出一種方向圖的旁瓣誤差模型用以評估性能，以尋找最優(yōu)參數(shù)設置. 下面以8×8元陣列為例進行說明.

首先定義水平視角和法線垂直視角兩個方向從歸一化零值開始搜索最小值和變化最大值，取這兩個值中較小的一個，作為該指向方向的主瓣邊界，定義垂直和水平方向的主瓣邊界所在矩形區(qū)域為主瓣區(qū)間. 然后確定最小迭代次數(shù)N及最大迭代次數(shù)M，8×8元陣列元陣列的最大迭代次數(shù)一般不超1 000. 第k次迭代的方向圖旁瓣誤差模型Ek定義為主瓣區(qū)間外的電平值與預期電平值(即預設虛擬干擾值)的差值. 則

Ek=Γk-Γ0.

(9)

如圖3(a)所示為旁瓣誤差隨迭代次數(shù)變化的曲線，此處沒有采用絕對值和對數(shù)是因為曲線與預設虛擬干擾值Γ0有關，如將虛擬干擾值設置為-20 dB，則曲線會整體小于零.

(a) 不同虛擬干擾下的旁瓣誤差

(b) 不同覆蓋區(qū)域下的旁瓣誤差

(c) 不同覆蓋區(qū)域的旁瓣誤差的導數(shù)圖3 利用方向圖旁瓣誤差模型確定迭代次數(shù)

針對實際情況，模型還可以更改覆蓋區(qū)間，如3(b)圖所示為四種不同覆蓋區(qū)域的旁瓣誤差變化曲線(所有方位：方位角0°～180°、俯仰角0°～90°；水平視角：方位角0°～180°、俯仰角0°；垂直視角：方位角90°、俯仰角0°～90°；綜合區(qū)域：俯仰角0°～5°及85°～90°時方位角為0°～180°、俯仰角6°～84°時方位角為主瓣邊界內(nèi)). 可以看出，所有方位和水平視角的曲線較為相似，垂直視角則誤差較大. 此處建議采用綜合區(qū)域.

為了確定迭代次數(shù)，先對上述曲線按照迭代次數(shù)作差，得到如圖3(c)所示的變化率曲線. 規(guī)定若滿足式(10)的條件，則認為曲線變化率已趨于平滑，此時迭代終止.

(10)

表2所示為設置不同最低門限時，迭代終止時的次數(shù).

表2 不同最低門限時的迭代終止值

由式(9)可以看出，此模型實際上是不考慮預設虛擬干擾值的絕對影響的，而只是針對某個固定虛擬干擾值時的性能波動. 在對預設虛擬干擾值進行設定時，需要對模型進行功能擴展.

2.2 實測數(shù)據(jù)分析

為檢驗算法對于實測數(shù)據(jù)的效果，利用上述方法對位于福建龍海的中程高頻地波雷達OSMAR071觀測數(shù)據(jù)進行分析. 高頻地波雷達OSMAR071是武漢大學電波傳播實驗室研制的陣列式OSMAR系列高頻地波雷達的定型產(chǎn)品，回波信號質(zhì)量較OSMAR工程樣機有較大的改善[16]. 該雷達是采用FMICW波形的寬波束地波雷達，波束范圍可達150°，掃頻帶寬為30 kHz，距離分辨率為5 km，接收陣列為8元非線性陣列，采用基于海洋回波的無源校準方法進行通道幅度和相位一致性校準后，利用MUSIC超分辨算法對海流進行測向，精度可達1.5°，該雷達每0.652 8×1 024 s可以得到一場多普勒譜數(shù)據(jù).

圖4為龍海雷達站2009年9月28日及2009年10月29日全天的實測數(shù)據(jù),每個時間點的數(shù)據(jù)為該時刻回波二階譜的頻率疊加，并對每個時間點的數(shù)據(jù)分別作二維CBF和二維NPS波束形成. 圖4(a)原始數(shù)據(jù)中4-6點時段存在點狀干擾，10-18點時段存在條狀干擾. 波束形成的主瓣指向陣列法線方向，工作頻率為7.88 MHz，主瓣指向陣列法線方向(方位角90°，俯仰角0°)；二維NPS算法設置旁瓣峰值為：-15 dB(俯仰角0°～5°)，-10 dB(俯仰角6°～90°)，迭代次數(shù)k=1 000次，每隔一個方位角、一個俯仰角取一個虛擬干擾，迭代的常數(shù)因子為η=0.1.

(a) 2009年9月28日全天的實測數(shù)據(jù)

(b) 2009年10月29日全天的實測數(shù)據(jù)圖4 龍海雷達站回波能量圖

比較結(jié)果可以發(fā)現(xiàn)，4-6點時段的點狀干擾數(shù)量明顯減少，10-18點時段的條狀干擾顏色變淺. 同樣的，圖4(b)也顯示出相同的性能. 這說明二維NPS能比二維CBF更有效地抑制點狀干擾，并能夠減輕條狀干擾的影響.

3 結(jié) 論

提出利用二維靜態(tài)方向圖NPS對高頻地波雷達回波信號進行數(shù)字波束形成，實現(xiàn)抑制電離層干擾的功能. 1)對于兩種陣型進行了仿真分析，發(fā)現(xiàn)對于8×8元陣列，二維NPS算法性能優(yōu)良，無論是在水平方向還是俯仰角較大的區(qū)域，都能夠抑制旁瓣電平，降低干擾噪聲對信號的影響；對于OSMAR071的小型陣列，需要仔細設置參數(shù)，特別是虛擬干擾功率大小，才能使綜合性能最優(yōu). 2)提出了方向圖旁瓣性能模型用來評估二維NPS算法迭代過程中的性能變化. 3)分析了兩天的實測數(shù)據(jù)，結(jié)果表明二維NPS算法能夠明顯減小電離層干擾的影響. 在實際應用中展示出較為穩(wěn)定可靠的性能.

由于實際電離層干擾的復雜性，抗干擾陣列的設計和算法的參數(shù)配置等需要進一步的研究. 今后還需要從以下問題開展進一步的工作：1)確定二維主瓣寬度的方法需要進行更多的比較嘗試. 2)方向圖性能模型的性能全面化，考慮抖動誤差和絕對誤差等多方面因素，便于迭代次數(shù)及虛擬干擾的選擇. 3)由于現(xiàn)有實測數(shù)據(jù)中電離層干擾都處于300～400 km的較遠區(qū)域，難以布設浮標獲取對比數(shù)據(jù). 同時海洋信息參數(shù)的獲取率很低，無法進行海洋參數(shù)的結(jié)果比較. 應設法尋找其它形式的電離層干擾進行實測數(shù)據(jù)分析.

[1] 吳海鵬, 焦培南, 凡俊梅. 高頻海洋回波譜電離層污染及實驗研究[J]. 電波科學學報, 2004, 19(5): 515-518.

WU HaiPeng, JIAO Peinan, FAN Junmei. Effects of ionosphere on HF sea echo Doppler spectra[J]. Chinese Journal of Radio Science, 2004, 19(5): 515-518. (in Chinese)

[2] LI W L. High-frequency over the horizon radar and ionospheric backscatter studies in China[J]. Radio Science, 1998,33(5):1445-1458.

[3] VOURAS P, FREBURGER B. Application of adaptive beamforming techniques to HF radar[C]∥ IEEE Radar Conference, 2008: 795-800.

[4] STOICA P, NEHORAI A. MUSIC, Maximum likelihood, and Cramer-Rao bound[J]. IEEE Trans on Acoustics, Speechand Signal Processing, 1989, 37(5): 720-741.

[5] 李 倫, 吳雄斌, 李 炎, 等. “鳳凰”臺風期間海面風、浪的高頻地波雷達OSMAR071遙感[J]. 遙感學報, 2012, 16(1): 154-165.

[6] 王永良, 丁前軍, 李榮鋒. 自適應陣列處理[M]. 北京: 清華大學出版社, 2009: 223-226.

[7] DOLPH C L. A Current distribution for broadside arrays which optimizes the relationship between beamwidth and side-lobe level[J]. Proc IRE, 1946, 34: 335-348.

[8] TSENG F I, CHENG D K. Optimum scannable planar arrays with an invariant sidelobe level[J]. proceedings of the IEEE, 1968, 56 (11): 1771-1778.

[9] BAKLANOV Y V. Chebyshev distribution of currents for a plane array of radiators[J]. Radio Eng Electron Phys, 1966, 11: 640-642.

[10] OLEN C A, COMPTON R T Jr. A numerical pattern synthesis algorithm for arrays[J]. IEEE Transactions on Antennas and Propagation, 1990, 38(10): 1666-1676.

[11] ER M H. Array pattern synthesis with a controlled mean-square sidelobel level[J]. IEEE Transactions on Signal Processing, 1992, 40(4): 977-981.

[12] SON S H, EOM S Y, JEON S I, et al. Automatic phase correction of phased array antennas by a genetic algorithm[J]. IEEE Antenas and Propagation Magazine, 2008, 56(8): 2751-2754.

[13] VILLEGAS F J. Parallel genetic-algorithm optimization of shaped beam coverage areas using planar 2-D phased arrays[J]. IEEE Trans on Antenna and Propagation, 2007, 55(6): 1745-1753.

[14] 尚 尚, 張 寧, 李 楊. 高頻地波雷達電離層雜波統(tǒng)計特性研究[J]. 電波科學學報, 2011, 26(3): 521-526.

SHANG Shang, ZHANG Ning, LI Yang. Ionospheric clutter statistical properties in HFSWR[J]. Chinese Journal of Radio Science,2011,26(3):521-526.(in Chinese)

[15] RAFAEL C, GONZALEZ. 數(shù)字圖像處理[M]. 北京: 電子工業(yè)出版社, 2002.

[16] WU X B, LI L, SHAO Y X, et al. Experimental determination of sinificant waveheight by OSMAR71: comparison with results from Buoy[J]. Journal of Wuhan University Nature Science, 2009, 14(6): 499-504.