雷達地形遮擋盲區的大氣折射修正算法

李冠運,劉松濤,徐華志

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

0 引言

地形遮蔽盲區的計算是雷達威力研究的重要內容。由于地形遮擋對雷達中低空探測產生的影響是基礎且廣泛的[1],研究雷達地形遮蔽盲區能夠為雷達選址、預警組網、低空突防等提供現實支撐,具有很好的研究應用價值。

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

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

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

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

1 大氣折射對雷達地形遮蔽盲區計算的影響

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

雷達地形遮蔽盲區是指在地形遮擋的作用下,雷達能量無法實現空間上的全面覆蓋,導致在雷達探測范圍內產生的不可探測空間域,如圖1所示。

圖1 地形遮擋導致的雷達探測盲區Fig.1 Radar detection blind spots due to terrain occlusion

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

在無線電領域,大氣折射通常指由于地球大氣層折射指數n在垂直梯度上的變化,造成無線電波傳播射線變得彎曲的現象。大氣折射計算出發點通常是Snell定律,以雷達為例,它表示為

nrcosθ=n0r0cosθ0=常數,

(1)

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

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

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

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

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

2.1 求解思路

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

圖4 遮擋關系示意圖Fig.4 Schematic diagram of shading relationship

2.2 主要問題及解決方法

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

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

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

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

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

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

(2)

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

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

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

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

根據等效地球半徑系數的計算公式:

(3)

得到k1～9≈1.21,k9～12≈1.05,對應的等效地球半徑分別為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,即可實現分級修正。

2) 具體算例

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

圖6 算例遮擋關系圖Fig.6 Arithmetic masking relationship diagram

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

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

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

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

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

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

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

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

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

2.3 計算結果的分析比較

使用2.2節中的算例,將未做大氣折射修正的算法記作f0(l),將使用單層等效地球法修正的算法記作fd1(l),將分層修正的算法記作fd2(l),分別計算地面投影距離lC取[16,400]時各算法計算結果,對比結果見圖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調整為800 m,以增大雷達初始仰角,再次計算f0(l),fd1(l)和fd2(l),對比結果見圖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)

設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,結果見圖13。

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

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

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

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

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

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

3.1 求解思路

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

(4)

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

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

圖14 平面直角坐標系構筑圖Fig.14 Image of the construction of a planar rectangular coordinate system

2) 求解O′的表達式

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

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

3) 求解h3

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

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

3.2 計算結果的分析比較

仍然使用2.2節中的算例,通過調節h2的取值來實現仰角的變化,評估其對fq(l)和fd(l)計算結果的影響(由于fq(l)無法顧及折射梯度的變化,因此這里選用fd1(l)做計算對比,此處選取標準大氣折射下的ae=8 490 km)。設θ0=0的方程算法為fq1(l),實際考慮θ0修正曲率半徑的方程算法為fq2(l),lC取[16～400],不同仰角的fd(l),fq1(l)和fq2(l)的計算結果對比見圖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θ在近零區間變化的鈍性,導致ρ的變化非常不明顯。因此,fd(l),fq1(l)和fq2(l)的計算結果差異很小。

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

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

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

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

4 結論

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

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

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

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