張 昊 王 璐 王文益 盧 丹 賈瓊瓊 吳仁彪

(中國(guó)民航大學(xué)天津市智能信號(hào)處理重點(diǎn)實(shí)驗(yàn)室，天津 300300)

1 引言

GNSS(Global Navigation Satellite System)信號(hào)到達(dá)用戶接收機(jī)時(shí)極其微弱,其信號(hào)強(qiáng)度在-160 dBW左右,極易受到各類干擾的影響。其中,欺騙式干擾與真實(shí)衛(wèi)星信號(hào)在時(shí)間域和頻率域都具有相同的格式,是一種故意的干擾。此外,由于欺騙式干擾的功率遠(yuǎn)低于接收機(jī)噪聲水平,傳統(tǒng)的干擾抑制方法不能很好地工作。因此,需要針對(duì)欺騙式干擾的特點(diǎn)利用新的方法對(duì)它進(jìn)行有效地抑制,保障GNSS接收機(jī)的定位、導(dǎo)航、授時(shí)正確[1-7]。

在現(xiàn)有的欺騙式干擾抑制方法中,多天線技術(shù)利用欺騙式干擾和真實(shí)衛(wèi)星信號(hào)空間特性的差異在空域上把欺騙式干擾和真實(shí)衛(wèi)星信號(hào)區(qū)分開,從而進(jìn)行欺騙式干擾的抑制。加拿大Calgary大學(xué)PLAN實(shí)驗(yàn)室的Saeed Daneshmand等人利用欺騙式干擾來自空間同一個(gè)方向的特點(diǎn),在空間能量最大的方向設(shè)置零陷來抑制欺騙式干擾,該方法不需要對(duì)陣列進(jìn)行校準(zhǔn),因此可以作為一個(gè)單獨(dú)的模塊安裝在接收機(jī)中[8]。Konovaltsev等人提出一種改進(jìn)的聯(lián)合欺騙檢測(cè)和姿態(tài)估計(jì)算法,通過波束形成技術(shù)在估計(jì)得到的欺騙式干擾來向上置零陷來抑制欺騙式干擾,并通過仿真和實(shí)采實(shí)驗(yàn)驗(yàn)證了欺騙式干擾抑制的性能[9]。多天線技術(shù)在欺騙式干擾抑制中已被證明是非常有效的方法,然而它的干擾抑制能力受自由度限制,而且多天線間通道誤差的存在將降低陣列處理性能。

合成陣列避免了傳統(tǒng)天線陣列的上述問題,其核心思想是把移動(dòng)天線接收到的時(shí)間采樣信號(hào)轉(zhuǎn)換成空間采樣信號(hào),通過單個(gè)天線運(yùn)動(dòng)形成合成陣列的多個(gè)虛擬陣元,并同時(shí)處理不同時(shí)間點(diǎn)不同空間位置上的采樣信號(hào)。合成陣列不存在通道一致性問題,同時(shí),天線運(yùn)動(dòng)可以形成更大的孔徑,使陣列的空間分辨率更高,因而來向估計(jì)和干擾抑制性能將更好。近年來,在衛(wèi)星導(dǎo)航領(lǐng)域應(yīng)用合成陣列的研究正逐漸開展。加拿大Calgary大學(xué)PLAN實(shí)驗(yàn)室的A.Broumandan等人基于軌跡精確可控的運(yùn)動(dòng)天線利用UCA-MUSIC算法對(duì)較大功率的壓制式干擾進(jìn)行了來向估計(jì)[10]。該實(shí)驗(yàn)室也提出了基于ESPRIT算法的AOA估計(jì)方法,并對(duì)真實(shí)的CDMA信號(hào)進(jìn)行了測(cè)試,在合成陣列的ESPRIT算法中,唯一的限制是陣元雙動(dòng)的平移不變性,此外,考慮到窄帶信號(hào)假設(shè)的影響,合成陣列的孔徑應(yīng)設(shè)定在幾米以內(nèi)[11]。該實(shí)驗(yàn)室還研究了GNSS信號(hào)的合成陣列穩(wěn)健波束形成算法,比較了多種波束形成方法的性能,但是,該方法需要已知GNSS信號(hào)的來向[12]。佛羅里達(dá)大學(xué)的Andrey Soloviev,俄亥俄大學(xué)的Frank van Graas以及空軍研究實(shí)驗(yàn)室的Mikel Miller研究了利用GPS天線運(yùn)動(dòng)合成大的天線孔徑產(chǎn)生窄波束的特點(diǎn)來進(jìn)行干擾消除[13]。

上述方法均是衛(wèi)星信號(hào)和本地碼相關(guān)之后再進(jìn)行處理,屬于基于相關(guān)域的方法,因此均需要改變接收機(jī)的結(jié)構(gòu),復(fù)雜度較高。此外,上述方法沒有針對(duì)欺騙式干擾進(jìn)行專門的研究。本文針對(duì)欺騙式干擾,在數(shù)據(jù)域直接處理,在不需要改變普通接收機(jī)結(jié)構(gòu)的前提下,提出一種利用合成陣列進(jìn)行GNSS欺騙式干擾抑制的方案。文章的結(jié)構(gòu)組織如下,第二部分推導(dǎo)了合成陣列接收數(shù)據(jù)模型,第三部分給出了數(shù)據(jù)域欺騙式干擾抑制的方法,第四部分為仿真實(shí)驗(yàn),第五部分為本文的總結(jié)。

2 數(shù)據(jù)模型

考慮一個(gè)軌跡精確可控的天線沿直線勻速運(yùn)動(dòng)如圖1所示。

圖1 合成陣列示意圖

圖1中x軸正方向表示單天線的運(yùn)動(dòng)方向,pant,1、pant,m表示天線在t1時(shí)刻和tm時(shí)刻的位置矢量,m表示第m個(gè)虛擬陣元,M表示虛擬陣元數(shù),d表示兩個(gè)相鄰虛擬陣元之間的距離,一般取半波長(zhǎng)。y軸為陣元運(yùn)動(dòng)方向的法線方向,psat,1和psat,m分別表示衛(wèi)星在t1時(shí)刻和tm時(shí)刻的位置矢量,ps表示欺騙式干擾的位置矢量,ρ1和ρm分別表示t1時(shí)刻和tm時(shí)刻衛(wèi)星和天線之間的距離,θ表示衛(wèi)星信號(hào)入射角度。

令t1時(shí)刻天線接收的衛(wèi)星信號(hào)為

s(t1)=Ad(t1-τ)c(t1-τ)exp(jωRF(t1-τ))

(1)

A表示衛(wèi)星信號(hào)的幅度,d(t)、c(t)分別代表衛(wèi)星信號(hào)的導(dǎo)航電文和擴(kuò)頻碼,ωRF表示射頻角頻率,則tm時(shí)刻天線接收的衛(wèi)星信號(hào)可寫作

s(tm)=Ad(tm-τ-Δτm)c(tm-τ-Δτm)

exp{jωRF(tm-τ-Δτm)}=

Ad(tm-τ-Δτm)c(tm-τ-Δτm)

exp{jωRF(tm-τ)+Δφm)}

(2)

其中,tm=t1+(m-1)ΔT,ΔT表示相鄰虛擬陣元的時(shí)間間隔。假設(shè)ΔT為C/A碼周期的整數(shù)倍(本文取1 ms),且假設(shè)導(dǎo)航電文不翻轉(zhuǎn),此時(shí)各陣元接收到的信號(hào)包絡(luò)不發(fā)生變化,因此,兩個(gè)虛擬陣元的接收信號(hào)只差一個(gè)相位為Δφm的指數(shù)項(xiàng),即式(2)可寫為

s(tm)=Ad(tm-τ-Δτm)c(tm-τ-Δτm)

exp{jωRF(tm-τ-Δτm)}=Ad(t1-τ)

c(t1-τ)exp{jωRF(tm-τ)}×exp(jΔφm)

(3)

令Δτm和Δφm表示衛(wèi)星和天線相對(duì)運(yùn)動(dòng)過程中,tm時(shí)刻與t1時(shí)刻天線接收到該顆衛(wèi)星的碼延遲之差和載波相位之差,Δτm和Δφm可以寫成如下形式

(4)

Δφm=-ωRFΔτm=-2π(ΔPsv-ΔPant)·ulos/λRF=

Δφsv,m+Δφant,m

(5)

因此,tm時(shí)刻天線接收的衛(wèi)星信號(hào)可重新寫作

s(tm)=s(t1)·exp(jΔφsv,m)·

exp(-jωd,sv(m-1)ΔT)

(6)

圖2利用高保真衛(wèi)星信號(hào)仿真器仿真了運(yùn)動(dòng)天線和靜止天線接收的衛(wèi)星信號(hào)多普勒頻率的變化情況,其中,運(yùn)動(dòng)天線為沿軌跡精確可控的軌道以20 m/s的速度作勻速直線運(yùn)動(dòng)。

圖3 僅由天線運(yùn)動(dòng)引起的多普勒頻率

由圖2可見,短時(shí)間內(nèi)當(dāng)天線固定不動(dòng)時(shí),多普勒頻率(此時(shí)僅由衛(wèi)星運(yùn)動(dòng)引起)呈線性變化,運(yùn)動(dòng)天線的多普勒頻率也呈線性變化且變化趨勢(shì)與靜態(tài)天線相同,兩者差值表示由天線運(yùn)動(dòng)引起的多普勒變化,且值為常數(shù)(如圖3所示),該值與天線的運(yùn)動(dòng)速度有關(guān)。因此,由天線運(yùn)動(dòng)引起的相位Δφant,m呈線性變化,和傳統(tǒng)陣列中陣列流型矢量的相位項(xiàng)類似。

3 數(shù)據(jù)域欺騙式干擾抑制

數(shù)據(jù)域欺騙式干擾抑制方案分為以下幾個(gè)模塊:合成陣列數(shù)據(jù)接收模塊、干擾來向估計(jì)模塊、天線權(quán)值計(jì)算模塊以及欺騙式干擾抑制模塊,如圖4所示。

圖4 數(shù)據(jù)域欺騙式干擾抑制算法原理框圖

3.1 合成陣列接收數(shù)據(jù)

合成陣列數(shù)據(jù)接收模塊由單個(gè)天線沿精確可控的固定軌跡勻速直線運(yùn)動(dòng)形成,天線在運(yùn)動(dòng)過程中接收衛(wèi)星信號(hào)和欺騙式干擾信號(hào),相鄰兩個(gè)虛擬陣元以ΔT為時(shí)間間隔接收采樣信號(hào)。在ΔT時(shí)間間隔內(nèi)的快拍數(shù)用字母N表示,在存在欺騙式干擾環(huán)境下,將單天線接收的信號(hào)經(jīng)下變頻和A/D采樣后的數(shù)字中頻信號(hào)寫成矢量形式,即可得到合成陣列的信號(hào)模型

(7)

3.2 干擾來向估計(jì)模塊

在得到合成陣列接收數(shù)據(jù)模型的基礎(chǔ)上,把合成陣列接收到的中頻采樣信號(hào)送入干擾來向估計(jì)模塊中分別利用MUSIC算法和RELAX算法進(jìn)行干擾來向估計(jì)。

3.2.1利用MUSIC算法估計(jì)干擾來向

多重信號(hào)分類(Multiple Signal Classification, MUSIC)是最為經(jīng)典的高分辨率來向估計(jì)算法,是一種基于矩陣特征分解的方法,通過在空域內(nèi)進(jìn)行一維譜峰搜索得到信源DOA的估計(jì)值。已知信號(hào)源的個(gè)數(shù),MUSIC算法的步驟歸納如下:

根據(jù)接收到信號(hào)的N個(gè)快拍數(shù)據(jù)計(jì)算信號(hào)協(xié)方差矩陣R的估計(jì)值

(8)

N表示ΔT時(shí)間間隔內(nèi)的快拍數(shù),代入合成陣列接收數(shù)據(jù)模型(7)可得

(9)

其中,Pa、Ps分別表示真實(shí)信號(hào)功率、欺騙式干擾功率;σ2表示噪聲功率;Lq表示第q個(gè)欺騙式干擾源轉(zhuǎn)發(fā)的PRN個(gè)數(shù)?？紤]到欺騙式干擾相對(duì)真實(shí)衛(wèi)星信號(hào)的碼延遲通常大于一個(gè)碼片,且它們之間的空間相關(guān)性較弱,因此上式的表達(dá)中略去了真實(shí)衛(wèi)星信號(hào)和欺騙式干擾相關(guān)后的結(jié)果。

(10)

本方案中,MUSIC算法需要已知欺騙式干擾源數(shù)目從而構(gòu)造噪聲子空間,此外,欺騙式干擾和真實(shí)衛(wèi)星信號(hào)的功率接近且均遠(yuǎn)低于噪聲的特點(diǎn)會(huì)使得MUSIC算法可能無法準(zhǔn)確獲得噪聲子空間而導(dǎo)致性能變差。因此MUSIC算法在當(dāng)前應(yīng)用背景下性能不夠理想。此時(shí)考慮利用不需要知道信號(hào)源數(shù)目以及對(duì)信噪比不敏感的RELAX算法進(jìn)行干擾來向估計(jì)。

3.2.2利用RELAX算法估計(jì)干擾來向

(11)

上述最小二乘方程很難直接求解,解決此問題的常用方法是松弛算法,如RELAX算法[14-16]。RELAX算法把復(fù)雜的多維優(yōu)化問題轉(zhuǎn)化為一系列一維優(yōu)化問題,在估計(jì)出一個(gè)信號(hào)之后,RELAX算法把該信號(hào)減去作為新的xk(n)來估計(jì)下一個(gè)信號(hào),在得到下一個(gè)信號(hào)的估計(jì)值后,RELAX算法會(huì)把新估計(jì)出的信號(hào)作為已知量減去重新估計(jì)第一個(gè)信號(hào),重復(fù)上述過程反復(fù)迭代直到所有的信號(hào)估計(jì)值都實(shí)際收斂,從而得到全部欺騙干擾源的來向

(12)

(13)

RELAX算法對(duì)信噪比的變化不敏感且能準(zhǔn)確分辨出兩個(gè)來向相近的信號(hào),然而算法過程需要多次迭代運(yùn)算,計(jì)算較復(fù)雜[14]。

3.3 天線權(quán)值計(jì)算和干擾抑制模塊

在權(quán)值計(jì)算模塊中,利用上述方法估計(jì)出所有欺騙干擾源的來向后,就可以構(gòu)造得到天線自適應(yīng)權(quán)矢量

wopt=(I-A(θ)AH(θ)/M)·g

(14)

(15)

把處理后的數(shù)據(jù)送入普通接收機(jī)中即可得到正確的導(dǎo)航結(jié)果。

4 仿真實(shí)驗(yàn)

實(shí)驗(yàn)中運(yùn)動(dòng)天線沿精確可控軌跡勻速直線運(yùn)動(dòng)形成包含10個(gè)虛擬陣元的合成陣列,相鄰虛擬陣元時(shí)間間隔ΔT=1 ms,陣元間距為GNSS信號(hào)波長(zhǎng)的一半d=λ/2。真實(shí)衛(wèi)星信號(hào)為PRN1、PRN3、PRN22、PRN25,入射角度分別為0°、10°、-10°、-40°,SNR=-20 dB。兩個(gè)欺騙式干擾源均由PRN1、PRN3、PRN22、PRN25的轉(zhuǎn)發(fā)信號(hào)組成。

為了說明干噪比對(duì)RELAX和MUSIC兩種算法性能的影響,在INR=-18 dB以及INR=-10 dB條件下,比較了欺騙式干擾在入射角度間隔分別為兩個(gè)波束寬度,一個(gè)波束寬度以及2/3個(gè)波束寬度條件下的來向估計(jì)性能。本次實(shí)驗(yàn)中固定其中一個(gè)欺騙式干擾源的來向?yàn)?0°,另一個(gè)欺騙式干擾源的來向隨著所設(shè)干擾源來向間隔的變化而變化。

表1 INR=-10 dB干擾來向估計(jì)值均方根誤差(度)

表2 INR=-18 dB干擾來向估計(jì)值均方根誤差(度)

由表1可以看出,當(dāng)干噪比INR=-10 dB時(shí),隨著干擾來向間隔逐漸減小,兩種算法的干擾來向估計(jì)性能均變差,且MUSIC算法在該干噪比下的干擾來向估計(jì)性能要比RELAX算法好;當(dāng)把欺騙式干擾的干噪比降為INR=-18 dB,如表2所示,隨著干擾來向間隔逐漸縮小,兩種算法的干擾來向估計(jì)性能均變差,然而在該干噪比下,RELAX算法的干擾來向估計(jì)性能要優(yōu)于MUSIC算法,其中MUSIC算法的來向估計(jì)性能隨著干噪比的降低而嚴(yán)重下降,RELAX算法的來向估計(jì)性能隨著干噪比的降低下降緩慢,體現(xiàn)出了RELAX算法對(duì)干噪比不敏感的特點(diǎn)。在公開的文獻(xiàn)中,欺騙式干擾的功率值一般設(shè)定為比真實(shí)衛(wèi)星信號(hào)大2 dB左右,從而保障欺騙式干擾能成功的欺騙地面接收機(jī),因此,可以看到在本章數(shù)據(jù)域欺騙式干擾抑制方案中,RELAX算法是一種更好的選擇。

下面給出上述兩種算法在干噪比為-18 dB條件下干擾抑制的性能,以干擾間隔為兩個(gè)波束寬度為例,天線方向圖如圖5所示。

圖5 天線方向圖

可以看到,采用MUSIC算法和RELAX算法抑制干擾時(shí),陣列天線在兩個(gè)欺騙式干擾來向上均形成了較深的零陷,保證能把欺騙式干擾抑制掉。在其他方向上的天線增益較為平坦,保證真實(shí)衛(wèi)星信號(hào)能以很小的損失通過。下面以PRN1為例,通過PRN1的二維搜索結(jié)果圖來給出本文算法下欺騙式干擾抑制前后的效果。

圖6(a)可以看出,欺騙式干擾抑制前,接收機(jī)二維搜索捕獲結(jié)果中存在三個(gè)明顯的相關(guān)峰,由于轉(zhuǎn)發(fā)式欺騙干擾功率高于真實(shí)衛(wèi)星信號(hào)功率,因此兩個(gè)較高的峰對(duì)應(yīng)欺騙式干擾相關(guān)峰,較低的峰對(duì)應(yīng)真實(shí)衛(wèi)星信號(hào)相關(guān)峰。利用本文所提合成陣列欺騙式干擾抑制算法進(jìn)行干擾抑制后的捕獲結(jié)果如圖6(b)、(c),可以看出欺騙式干擾均已被成功抑制,接收機(jī)此時(shí)可以正常捕獲到真實(shí)衛(wèi)星信號(hào)。

圖6 PRN1捕獲結(jié)果對(duì)比

5 結(jié)論

本文基于合成陣列的場(chǎng)景,提出了一種欺騙式干擾抑制方案,首先分析了存在欺騙式干擾的環(huán)境下合成陣列的接收信號(hào)模型,文中利用RELAX算法和MUSIC算法估計(jì)了多個(gè)欺騙式干擾來向,并在該方案中比較了這兩種方法的來向估計(jì)性能,在干擾來向估計(jì)之后通過構(gòu)造合成陣列自適應(yīng)權(quán)向量來抑制欺騙式干擾并從接收機(jī)捕獲結(jié)果比較了這兩種算法的干擾抑制性能。上文所述方法針對(duì)天線運(yùn)動(dòng)軌跡精確可控,即合成陣列的陣列流型已知的情況,不考慮天線運(yùn)動(dòng)誤差的影響,而在實(shí)際應(yīng)用中,由陣元運(yùn)動(dòng)引起的位置誤差是不可避免的,因此為了體現(xiàn)合成陣列的優(yōu)點(diǎn)需要對(duì)陣元運(yùn)動(dòng)引起的位置誤差進(jìn)行補(bǔ)償,關(guān)于合成陣列運(yùn)動(dòng)誤差補(bǔ)償?shù)膯栴}將在下一步工作中進(jìn)行。

[1] 劉春保. 衛(wèi)星導(dǎo)航系統(tǒng)脆弱性評(píng)估與對(duì)策[J]. 衛(wèi)星應(yīng)用, 2015(4): 49-54.

Liu Chunbao. Vulnerability Assessment and Countermeasure of Satellite Navigation System[J]. Satellite Applications, 2015(4): 49-54. (in Chinese)

[2] 辛潔, 趙偉, 張之學(xué),等. 衛(wèi)星導(dǎo)航系統(tǒng)發(fā)展及其軍事應(yīng)用特點(diǎn)分析[J]. 導(dǎo)航定位學(xué)報(bào),2015(4):38- 43.

Xin Jie, Zhao Wei, Zhang Zhixue, et al. Development of GNSS and Characteristic Analysis of Military Application[J]. Journal of Navigation and Positioning, 2015(4): 38- 43. (in Chinese)

[3] 王群. 四大衛(wèi)星導(dǎo)航系統(tǒng)及其特點(diǎn)[J]. 國(guó)防科技, 2015(1): 40- 43.

Wang Qun. Four satellite navigation systems and their characteristics[J]. National Defense Science and Technology, 2015(1): 40- 43. (in Chinese)

[4] 成躍進(jìn). 現(xiàn)代衛(wèi)星導(dǎo)航定位系統(tǒng)發(fā)展介紹[J]. 空間電子技術(shù), 2015(1): 17-25.

Cheng Yuejin. The Development Intro of Satellite Navigation and Positioning System[J]. Space Electronic Technology, 2015(1): 17-25. (in Chinese)

[5] 盧丹. 穩(wěn)健的全球?qū)Ш叫l(wèi)星系統(tǒng)抗干擾技術(shù)研究[D]. 西安: 西安電子科技大學(xué), 2013.

Lu Dan. Robust Anti-jamming Techniques for Global Navigation Satellite System[D]. Xi’an: XiDian Univerisity, 2013.(in Chinese)

[6] Shen Jun, Ma Jiaqing, Yang Jian. Updated on Development of BeiDou Navigation Satellite System[C]∥Proceedings of the 28th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+2015), Tampa, Florida, USA, 2015: 1264-1290.

[7] Connor L. Brashar. An Analysis of Near-band Harmful Interference on Civilian GPS Receivers[C]∥Proceedings of the 29th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+2016), Portland, Oregon, 2016: 3123-3132.

[8] Daneshmand S, Jafarnia-Jahromi A, Broumandan A, et al. A Low-Complexity GPS Anti-Spoofing Method Using a Multi-Antenna Array[C]∥Proceedings of the 25th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2012), Nashville, TN, USA, 2012: 1233-1243.

[9] Konovaltsev A, Cuntz M, Haettich C, et al. Autonomous Spoofing Detection and Mitigation in a GNSS Receiver with an Adaptive Antenna Array[C]∥Proceedings of the 26th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+2013), Nashville, TN, USA, 2013: 2937-2948.

[10] Broumandan A, Lin Tao, Moghaddam A, et al. Direction of Arrival Estimation of GNSS Signals Based on Synthetic Antenna Array[C]∥Proceedings of the 20th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2007), Fort Worth, TX, 2007: 728-738.

[11] Ali Broumandan, John Nielsen, Gérard Lachapelle. Practical Results of High Resolution AOA Estimation by the Synthetic Array[C]∥Proceeding of IEEE Vehicular Technology Conference’08 Fall, Calgary, Alberta, Canada, Sep.21-24, 2008.

[12] Lin Tao, Ali Broumandan, John Nielsen, et al. Robust Beamforming for GNSS Synthetic Antenna Arrays[C]∥Proceedings of the 22th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2009), Savannah, GA, 2009:387- 401.

[13] Andrey Soloviev, Frank van Graas, Sanjeev Gunawardena, et al. Synthetic Aperture GPS Signal Processing: Concept and Feasibility Demonstration[C]∥Proceedings of the 2009 International Technical Meeting of the Institute of Navigation, Anaheim, CA, 2009: 851- 863.

[14] Li Jian, Petre Stoica. Efficient mixed-spectrum estimation with applications to target feature extraction[J]. IEEE Transactions on Signal Processing, 1996, SP- 44(2): 281-295.

[15] Li Jian, Zheng Dunmin, Petre Stoica. Angle and waveform estimation via RELAX[J]. IEEE Transactions on Aerospace and Electronic Systems, 1997, 33(3): 1077-1087.

[16] Liu Zhengshe, Li Jian. Implementation of the RELAX algorithm[J]. IEEE Transactions on Aerospace and Electronic Systems, 1998, 34(2): 657- 664.