李冠運(yùn),劉松濤,徐華志

(1.中國人民解放軍92011部隊,上海 202150;2.海軍大連艦艇學(xué)院信息系統(tǒng)系,遼寧 大連 116018)

0 引言

地形遮蔽盲區(qū)的計算是雷達(dá)威力研究的重要內(nèi)容。由于地形遮擋對雷達(dá)中低空探測產(chǎn)生的影響是基礎(chǔ)且廣泛的[1],研究雷達(dá)地形遮蔽盲區(qū)能夠為雷達(dá)選址、預(yù)警組網(wǎng)、低空突防等提供現(xiàn)實(shí)支撐,具有很好的研究應(yīng)用價值。

目前,地理信息系統(tǒng)(GIS)技術(shù)的發(fā)展與應(yīng)用為量化評估地形遮蔽對雷達(dá)探測造成的影響提供了研究條件,但目前的相關(guān)結(jié)合性應(yīng)用研究多依賴GIS技術(shù)的通視分析功能,主要有基于視線的方法、基于并行計算系統(tǒng)的可視性分析模型、利用鄰近視點(diǎn)之間的內(nèi)在一致性改進(jìn)并行算法、基于斜率比較的可視性算法等[2]。然而,地形遮蔽條件下的雷達(dá)探測能力分析并不完全等同于通視分析,若直接套用相關(guān)算法將導(dǎo)致以下兩個問題。

一是計算結(jié)果不精確。通視分析相關(guān)算法多聚焦計算效率和準(zhǔn)確性的提升,往往不考慮大氣折射誤差的修正。然而對于雷達(dá)探測來說,大氣折射將導(dǎo)致雷達(dá)波束彎曲,顯著影響遮蔽盲區(qū)的計算,若直接套用通視算法而不修正大氣折射,將導(dǎo)致計算出現(xiàn)嚴(yán)重誤差。但是對相關(guān)誤差的修正研究并不多見。

二是計算產(chǎn)品的適用性不強(qiáng)。通視分析多聚焦于觀察點(diǎn)與地表或某高度層的可視分析,形成的產(chǎn)品多為可視域。文獻(xiàn)[3]即是在通視分析的基礎(chǔ)上,求解地形遮擋條件下的地表電磁波覆蓋范圍,并運(yùn)用等效地球法對結(jié)果進(jìn)行了修正。

不同于以上研究,本文所聚焦的是地表以上的全空域分析,目的是對雷達(dá)威力邊界進(jìn)行定位,計算產(chǎn)品更具有普適性和實(shí)時性。本文提出運(yùn)用等效地球半徑法和二元曲率方程法來計算雷達(dá)地形遮蔽時的大氣折射誤差,可以有效地提升雷達(dá)地形遮蔽盲區(qū)的建模精度和計算產(chǎn)品的適用性。

1 大氣折射對雷達(dá)地形遮蔽盲區(qū)計算的影響

在雷達(dá)探測領(lǐng)域,針對大氣折射導(dǎo)致的測角、測距、測速等誤差修正[4-5],已經(jīng)形成了較為成熟的方法。文獻(xiàn)[6]對射線描跡法、線性分層法和等效地球半徑法等誤差修正方法在多種場景下的具體應(yīng)用進(jìn)行了詳細(xì)說明。文獻(xiàn)[7]提出了利用差分方法求解任意大氣層中的射線規(guī)范方程,并對單脈沖雷達(dá)和三基地雷達(dá)進(jìn)行了實(shí)用的誤差修正。文獻(xiàn)[8]對各類算法的使用范圍及優(yōu)缺點(diǎn)進(jìn)行了分析,并給出了雷達(dá)系統(tǒng)大氣折射誤差修正技術(shù)在今后的研究方向。文獻(xiàn)[9]為解決誤差修正的精確性與實(shí)時性矛盾提供了思路。本文借鑒以上解算方法,研究雷達(dá)地形遮擋盲區(qū)受大氣折射的影響情況。

雷達(dá)地形遮蔽盲區(qū)是指在地形遮擋的作用下,雷達(dá)能量無法實(shí)現(xiàn)空間上的全面覆蓋,導(dǎo)致在雷達(dá)探測范圍內(nèi)產(chǎn)生的不可探測空間域,如圖1所示。

圖1 地形遮擋導(dǎo)致的雷達(dá)探測盲區(qū)Fig.1 Radar detection blind spots due to terrain occlusion

若不考慮Fresnel區(qū)影響,地物遮擋導(dǎo)致雷達(dá)威力范圍在三維空間中存在邊界,邊界以下為地形遮蔽盲區(qū)。計算雷達(dá)遮蔽盲區(qū)的本質(zhì)即是求解該邊界,以角度為切分,該邊界應(yīng)是以雷達(dá)距離l為變量的高度函數(shù),記為f(l)。

在無線電領(lǐng)域,大氣折射通常指由于地球大氣層折射指數(shù)n在垂直梯度上的變化,造成無線電波傳播射線變得彎曲的現(xiàn)象。大氣折射計算出發(fā)點(diǎn)通常是Snell定律,以雷達(dá)為例,它表示為

nrcosθ=n0r0cosθ0=常數(shù),

(1)

式(1)中,θ0和θ分別為射線初始仰角和射線仰角,n和n0分別為地面和空中折射指數(shù),r和r0分別為雷達(dá)和目標(biāo)到地心距離。定義e0為射線上某點(diǎn)切線指向角對弧長的變化率,曲率半徑ρ=1/e0,根據(jù)式(1)有ρ=1/e0=-n/[cosθ(dn/dh)],考慮到n≈1,可見折射指數(shù)垂直梯度dn/dh越大,射線曲率半徑ρ就越小,射線彎曲就越明顯[10]。大氣折射對計算雷達(dá)地形遮擋盲區(qū)造成的影響如圖2所示。

圖2 大氣折射對雷達(dá)地形遮擋盲區(qū)造成的影響Fig.2 The effect of atmospheric refraction on the blind spot caused by radar terrain occlusion

圖2中,雷達(dá)A波束被山B遮擋,形成了威力邊界f(l),由圖中不難看出,受大氣折射影響,f(l)向地心方向彎曲,導(dǎo)致雷達(dá)探測盲區(qū)高度降低。該誤差會隨著雷達(dá)作用距離l增大而快速提升,使求解結(jié)果嚴(yán)重失真,如圖3所示。因此,對大氣折射的修正是求解雷達(dá)地形遮擋盲區(qū)關(guān)鍵環(huán)節(jié)。

圖3 雷達(dá)地形遮擋盲區(qū)高度隨雷達(dá)作用距離變化的示意圖Fig.3 Schematic of the variation of radar terrain occlusion blind spot height with radar action distance

2 基于等效地球半徑的大氣折射修正算法

2.1 求解思路

等效地球半徑法是在一定條件下,將實(shí)際地球半徑以某種等效地球半徑替代,從而使大氣等效為均勻,實(shí)現(xiàn)化曲為直,簡化計算目的,在工程上得到廣泛引用。關(guān)于等效地球半徑法推導(dǎo)過程,文獻(xiàn)[10-11]均有較為細(xì)致闡明,在此不做贅述。等效地球半徑法修正大氣折射是在雷達(dá)威力邊界計算中引入了基于等效地球半徑的電磁波傳播模型來實(shí)現(xiàn),下面以實(shí)例簡要介紹該方法的應(yīng)用。

圖4 遮擋關(guān)系示意圖Fig.4 Schematic diagram of shading relationship

2.2 主要問題及解決方法

2.2.1初始仰角θ0可能引起的誤差

關(guān)于等效地球半徑法的概念定義和局限性問題,文獻(xiàn)[11-12,14]做了討論,但對于起始仰角θ0≈0是否為該方法的充分必要條件,給出了不同意見。文獻(xiàn)[11]以等效地球中相對曲率必須保持不變?yōu)橐罁?jù),得出了θ0≈0是等效地區(qū)半徑法充分必要條件;文獻(xiàn)[12]從大氣折射指數(shù)梯度近似角度出發(fā),認(rèn)為θ0取值不影響等效地球半徑法的適用;文獻(xiàn)[14]則認(rèn)為該問題主要取決于等效地球半徑法定義方法,若以相對曲率等效法來看θ0≈0是必須的,但若以折射梯度效應(yīng)等效法來推導(dǎo),θ0取值與方法的使用基本無關(guān)。本文通過數(shù)據(jù)計算的方式檢驗證明,運(yùn)用等效地球半徑法計算雷達(dá)地形遮蔽盲區(qū)這一工程應(yīng)用中,θ0≈0并不是等效地球法的限制條件。

2.2.2大氣折射梯度的非線性變化可能引起的誤差

等效地球半徑法的定義是基于均勻球面分層且折射指數(shù)垂直梯度不隨高度變化(dn/dh為常數(shù))的大氣環(huán)境[11]。根據(jù)文獻(xiàn)[10],在近地0～1 km高度范圍內(nèi),dn/dh呈線性分布,4/3等效地球半徑法精度較高,但高于1 km后dn/dh為負(fù)指數(shù)分布,導(dǎo)致等效地球半徑法出現(xiàn)誤差,而防空雷達(dá)地形遮擋盲區(qū)高度往往超過了這一數(shù)值。針對以上不足,在文獻(xiàn)[13]的基礎(chǔ)上,本文提出一種多層等效地球半徑法,實(shí)現(xiàn)對誤差的進(jìn)一步修正。

1) 確定多層等效地球半徑

根據(jù)文獻(xiàn)[10],地面至60 km高度范圍中,較為科學(xué)的大氣折射描述方法應(yīng)為分段模型,即海拔高度h處的折射率N(h)為

(2)

式(2)中,N0為地面折射率,h0為地面海拔,ΔN1是0～1 km折射變化率(標(biāo)準(zhǔn)大氣下是39 N/km),N1是海拔1 km折射率,Ca1為地面上空1～9 km指數(shù)衰減系數(shù),N9是海拔9 km折射率(該值一般穩(wěn)定在105 N),Ca9為地面上空9～60 km指數(shù)衰減系數(shù)(每千米年平均衰減系數(shù)為0.143 2)。

因大氣環(huán)境的差異,N0,ΔN1,Ca1在不同地區(qū)、不同時段取值均有不同,在計算時可參考各地區(qū)統(tǒng)計數(shù)據(jù)。在此以北京地區(qū)7月份平均統(tǒng)計數(shù)據(jù)為例,取N0=369 N,ΔN1=51 N/km,Ca1=0.14,繪制1～9 km和9～20 km負(fù)指數(shù)函數(shù)圖像,見圖5。

圖5 北京地區(qū)7月份1～20 km大氣折射率變化情況圖Fig.5 Changes in atmospheric refractive index in Beijing from January to July

圖5中,N(h)在1～20 km區(qū)間函數(shù)形態(tài)仍保留較強(qiáng)線性特征,若將其近似為線性函數(shù),可實(shí)現(xiàn)在高空域使用等效地球法的目的,大大簡化復(fù)雜運(yùn)算而付出較小的精度代價。根據(jù)計算,1～9 km平均ΔN=27 N/km;9～20 km平均ΔN≈7.5 N/km。

根據(jù)等效地球半徑系數(shù)的計算公式:

(3)

得到k1～9≈1.21,k9～12≈1.05,對應(yīng)的等效地球半徑分別為a1=7 708 km,a2=6 689 km。通過在0～1 km,1～9 km,9～20 km空域分別使ae等于8 490 km,7 708 km和6 689 km,即可實(shí)現(xiàn)分級修正。

2) 具體算例

假設(shè)雷達(dá)A架設(shè)高度h1=72 m,在某方向上受到山B遮擋,山B與雷達(dá)A相距為L2=15.14 km,山B高度h2=222 m,見圖6,求解其400 km外的C點(diǎn)地形遮蔽盲區(qū)高度h3。

圖6 算例遮擋關(guān)系圖Fig.6 Arithmetic masking relationship diagram

首先取ae=8 490 km做等效地球變換,則有圖7所示遮擋關(guān)系。

圖7 算例第一層遮擋關(guān)系圖Fig.7 Image of the first layer of occlusion relations of the algorithm

假定存在G點(diǎn)海拔為1 000 m,根據(jù)幾何換算關(guān)系,可反推出G點(diǎn)與雷達(dá)A的地表投影距離為70.5 km,G點(diǎn)與C點(diǎn)地表距離為329.5 km,h3=13 118 m。以G點(diǎn)為起點(diǎn),ae=7 708 km,做第二層等效地球轉(zhuǎn)化,在此最困難的是確定轉(zhuǎn)換后的幾何關(guān)系,解決方法是參考第一次轉(zhuǎn)化,幾何關(guān)系的確定是依賴山B的“錨點(diǎn)”作用,因此在第二層轉(zhuǎn)化前,要先找到G點(diǎn)的“錨點(diǎn)”。假設(shè)距離G點(diǎn)地表距離100 km存在高度為h4的山H,其對G點(diǎn)的遮擋剛好在C點(diǎn)達(dá)到13 118 m,根據(jù)幾何關(guān)系,h4=3 322 m,則第二層等效地球轉(zhuǎn)化可基于山H對G點(diǎn)的遮擋,見圖8。

圖8 算例第二層遮擋關(guān)系圖Fig.8 Image of the second layer of occlusion relations of the algorithm

進(jìn)一步假定存在J點(diǎn)海拔為9 000 m,根據(jù)幾何換算關(guān)系,可反推出J點(diǎn)與G點(diǎn)的地表投影距離為245.02 km,J與C點(diǎn)地表距離為84.48 km,h3=13 573 m。以J點(diǎn)為起點(diǎn),ae=6 689 km,做第三層等效地球轉(zhuǎn)化前假設(shè)距J點(diǎn)50 km的山K(高度h5=11 594 m)是本次轉(zhuǎn)化的“錨點(diǎn)”,則有以下遮擋關(guān)系,見圖9。

圖9 算例第三層等效地球變換前的遮擋關(guān)系圖Fig.9 Image of occlusion relations in front of the equivalent earth in the third layer of the algorithm

最后以山K遮擋J點(diǎn)為基礎(chǔ),做第三層等效地球轉(zhuǎn)化,見圖10,求得h3=13 601 m。

圖10 算例第三層等效地球變換后的遮擋關(guān)系圖Fig.10 Image of occlusion relations after the equivalent earth in the third layer of the algorithm

從修正結(jié)果看,h3的取值從13 118 m修正為13 601 m,代表了大氣折射的弱化,符合高空域大氣折射降低的實(shí)際。

2.3 計算結(jié)果的分析比較

使用2.2節(jié)中的算例,將未做大氣折射修正的算法記作f0(l),將使用單層等效地球法修正的算法記作fd1(l),將分層修正的算法記作fd2(l),分別計算地面投影距離lC取[16,400]時各算法計算結(jié)果,對比結(jié)果見圖11。

圖11 f0(l),fd1(l),fd2(l)隨lC變化趨勢對比圖(h2=222 m)Fig.11 Trends comparative chart,f0(l),fd1(l),fd2(l) with lC(h2=222 m)

為便于分析比對,將算例中的山B高度h2調(diào)整為800 m,以增大雷達(dá)初始仰角,再次計算f0(l),fd1(l)和fd2(l),對比結(jié)果見圖12。

圖12 f0(l),fd1(l),fd2(l)隨lC變化趨勢對比圖(h2=800 m)Fig.12 Trends comparative chart,f0(l),fd1(l),fd2(l) with lC(h2=800 m)

設(shè)t1(lC)=100%·[f0(lC)-fd1(lC)]/f0(lC),t2(lC)=100%·[f0(lC)-fd2(lC)]/f0(lC),其值可表示兩種算法與f0(l)的比例差值。計算h2=222 m和h2=800 m時的t1,t2,結(jié)果見圖13。

圖13 t1,t2隨lC變化趨勢對比圖Fig.13 Trends comparative chart,t1,t2 with lC

分析圖11-圖13,不難看出:

1) 隨著lC距離的增大和電磁波所在空域的提升,fd1(l)和fd2(l)之間無論是絕對差值或是比例差值,均呈現(xiàn)不斷擴(kuò)大趨勢。計算結(jié)果反映了大氣折射現(xiàn)象在高空域減弱這一事實(shí)。

2)fd1(l)和fd2(l)的差值與電磁波在高空域的傳播距離相關(guān),電磁波越早進(jìn)入1 km以上空域,兩者差值越早開始擴(kuò)大,但隨著電磁波在高空域長時間運(yùn)行,該差值有逐步趨穩(wěn)態(tài)勢。對于防空雷達(dá)關(guān)注的20 km以下空域,兩者的差值比較有限,若對精度要求不高,可以直接使用fd1(l)進(jìn)行計算。

3) 從t1,t2來看,等效地球法對高仰角的修正效果要明顯弱于低仰角的情況,這與Snell定律表述相符,即仰角越高大氣折射效應(yīng)越弱。這也反映了等效地球法對仰角變化具有敏感性。

3 基于二元曲率方程的大氣折射修正算法

3.1 求解思路

根據(jù)文獻(xiàn)[14],大氣折射導(dǎo)致雷達(dá)波束以曲線形式傳播,射線曲率e0定義為曲線上某點(diǎn)切線指向角對弧長的變化率。衡量射線的彎曲程度可用曲率半徑ρ=1/e0來表示,曲率半徑越大,代表彎曲程度越小,反之越大,其表達(dá)式為

(4)

式(4)中,n0為地面大氣折射指數(shù),a為地球半徑約6 370 km,θ0為射線的初始仰角,n為大氣高度h處的大氣折射指數(shù),dn/dh為折射指數(shù)垂直梯度。由此可見,射線的彎曲是由于存在dn/dh而產(chǎn)生的折射效應(yīng),由于n≈1,h?a,故當(dāng)dn/dh為恒定值情況下,射線曲率取決于起始角θ0,則曲率半徑ρ為恒定值,若將曲線無限延長,射線最終將首尾相接并構(gòu)成標(biāo)準(zhǔn)圓。若雷達(dá)A被山B遮擋,則該曲率圓必過A,B兩點(diǎn),在ρ已知的情況下圓的方程是可求的,繼而可修正大氣折射。具體解法為:

1) 以地球圓心O為原點(diǎn),構(gòu)建平面直角坐標(biāo),見圖14。假定存在雷達(dá)A高度為h1,山B高度為h2,山B遮擋雷達(dá)A導(dǎo)致在C點(diǎn)形成高度為h3的雷達(dá)盲區(qū),地球半徑為a,AB之間地表距離為lB,對應(yīng)地心角為α,AC對應(yīng)地表距離為lC,對應(yīng)地心角為α+β,大氣折射形成的曲率圓心為O′,曲率半徑為ρ。則有地球方程表達(dá)式為x2+y2=a2,A點(diǎn)坐標(biāo)為(0,a+h1),B點(diǎn)坐標(biāo)為(asinα,acosα)。h1,h2,lB,lC,a均為已知量,α,β可由已知量推算。

圖14 平面直角坐標(biāo)系構(gòu)筑圖Fig.14 Image of the construction of a planar rectangular coordinate system

2) 求解O′的表達(dá)式

首先確定曲率半徑ρ。標(biāo)準(zhǔn)大氣折射情況下ρS=4a[6],則ρ=ρS/cosθ0。由于ρ?lB,得θ0≈∠BAO-π/2,△ABO可解,故θ0,ρ可解。

設(shè)O′表達(dá)式為(x-b)2+(y-c)2=ρ2,將A、B坐標(biāo)代入,解二元二次方程,排除掉c>0的根項即得到O′的方程表達(dá)式。

3) 求解h3

直線OC的表達(dá)式為y=x·cot(α+β),將OC表達(dá)式與O′的方程表達(dá)式聯(lián)立,再次解二元二次方程,排除掉小于0的根,即可求得C點(diǎn)的坐標(biāo),繼而求得h3。文中將通過曲率方程法修正大氣折射的算法記為fq(l)。

可以看到,曲率方程法具有清晰的物理概念,計算過程簡單,能夠考慮θ0造成的誤差。缺點(diǎn)是解二元二次方程算力代價較大,且dn/dh必須為恒定值,高空大氣折射的非線性變化會帶來一定誤差。

3.2 計算結(jié)果的分析比較

仍然使用2.2節(jié)中的算例,通過調(diào)節(jié)h2的取值來實(shí)現(xiàn)仰角的變化,評估其對fq(l)和fd(l)計算結(jié)果的影響(由于fq(l)無法顧及折射梯度的變化,因此這里選用fd1(l)做計算對比,此處選取標(biāo)準(zhǔn)大氣折射下的ae=8 490 km)。設(shè)θ0=0的方程算法為fq1(l),實(shí)際考慮θ0修正曲率半徑的方程算法為fq2(l),lC取[16～400],不同仰角的fd(l),fq1(l)和fq2(l)的計算結(jié)果對比見圖15和圖16。

圖15 fd(l),fq1(l),fq2(l)取值計算圖(h2=2 000 m,θ=7.188°,ρ=25 681.84 km)Fig.15 Value calculation chart of fd(l),fq1(l),fq2(l)(h2=2 000 m,θ=7.188°,ρ=25 681.84 km)

圖16 fd(l),fq1(l),fq2(l)取值計算圖(h2=100 000 m,θ=81.25°,ρ=167 495 km)Fig.16 Value calculation chart of fd(l),fq1(l),fq2(l)(h2=100 000 m,θ=81.25°,ρ=167 495 km)

由圖15可知,由于cosθ在近零區(qū)間變化的鈍性,導(dǎo)致ρ的變化非常不明顯。因此,fd(l),fq1(l)和fq2(l)的計算結(jié)果差異很小。

當(dāng)θ0增大到一定程度后,fq1(l)與fq2(l)的計算差值逐漸明顯,而fd(l)的變化趨勢總體與fq2(l)相同,如圖16所示。

通過圖15和圖16不難得出,在起始仰角如此大幅變化的情況下,fq(l),fd(l)的計算結(jié)果呈現(xiàn)出高度擬合性,這至少能夠說明兩點(diǎn):

1) 由于fq(l)和fd(l)采用了完全不同求解思路,而最終計算結(jié)果的擬合性證實(shí)了兩種方法均為可信算法。

2) 由于fq(l)顧及了起始仰角對計算結(jié)果的影響,fq(l)與fd(l)的高度重合結(jié)合2.3節(jié)中fd(l)的敏感性結(jié)論,表明了θ0≈0并不是等效地球算法的限制條件。

4 結(jié)論

本文首次系統(tǒng)研究了雷達(dá)地形遮蔽盲區(qū)計算中大氣折射的修正方法,并針對大氣折射梯度非線性變化帶來的誤差,提出了借用“錨點(diǎn)”的概念實(shí)現(xiàn)分層等效地球的方法;同時,針對起始仰角θ0可能帶來的誤差,提出了構(gòu)建曲率方程的方法評估仰角變化對計算結(jié)果的影響。最后從多維度檢驗和驗證實(shí)現(xiàn)了大氣折射修正的等效地球半徑法和曲率方程法的互洽,結(jié)果表明:

1) 對于5 km以下的中低空大氣,單層等效地球半徑法和二元曲率方程法精確度基本能夠滿足防空雷達(dá)地形遮蔽盲區(qū)計算的需求;當(dāng)空域提升至5 km以上時,若對計算精度有較高要求,則需選用多層等效地球半徑法進(jìn)行修正。

2) 經(jīng)與二元曲率方程法對比檢驗,證明等效地球半徑法考慮了起始仰角對計算結(jié)果造成的擾動,可適用于起始角θ0≠π/2的計算條件。

3) 若以算力代價來看,等效地球半徑法不涉及高次元計算,計算速度較快。若以建模便捷性來看,二元曲率方程法物理概念明晰,不涉及對地球的變換,更適合GIS技術(shù)背景下的模型構(gòu)造和可視化呈現(xiàn)。