快速算法實(shí)現(xiàn)電離層短波射線追蹤

孫 方， 康士峰， 趙振維， 王紅光

(中國(guó)電波傳播研究所青島分所，山東 青島 266107)

0 引言

目前，在計(jì)算電離層短波射線追蹤方面，求解基于費(fèi)爾馬原理得到的射線微分方程，是計(jì)算射線軌跡的比較理想和普遍的數(shù)值方法。然而，在射線方程的求解過程中，提高精度與減少運(yùn)算時(shí)間成為一對(duì)矛盾：如果將群路徑步長(zhǎng)設(shè)的太大，將引起不可忽略的誤差，反之，計(jì)算時(shí)間將會(huì)成倍的增加。所以，人們研究了快速的射線追蹤算法，達(dá)到了精度降低很小的前提下，大大提高計(jì)算速度的目的。目前主要從如下兩個(gè)方面實(shí)現(xiàn)快速追蹤：①在電離層以下采用直線近似傳播的方式[1]；②采用變步長(zhǎng)技術(shù)，根據(jù)電離層等離子體頻率梯度的大小自適應(yīng)改變?nèi)郝窂讲介L(zhǎng)。結(jié)合這兩方面，詳細(xì)給出了快速算法的計(jì)算過程，并引入亞大模型和附加擾動(dòng)影響的準(zhǔn)拋物模型作為電離層剖面對(duì)計(jì)算結(jié)果進(jìn)行了比較，體現(xiàn)了快速算法實(shí)現(xiàn)電離層短波射線追蹤的優(yōu)越性。

1 快速算法

快速算法的主要思想是將射線路徑按傳播順序分成三段（如圖 1所示），分別是發(fā)射點(diǎn)到電離層底高（TA）段、電離層段（AB）以及電離層底高到地面（BG）段。

在TA和BG段，射線在電離層外，將射線看作直線傳播，利用幾何關(guān)系直接計(jì)算；而在AB段，則利用球坐標(biāo)下的射線方程進(jìn)行數(shù)值計(jì)算，并采用變步長(zhǎng)技術(shù)，使算法達(dá)到準(zhǔn)確且計(jì)算速度快的目的。下面將詳細(xì)分析各段的計(jì)算方法。

1.1 發(fā)射點(diǎn)到電離層底高段軌跡方程

由于電離層段使用的是球坐標(biāo)下的射線方程，為了方便銜接計(jì)算，這里的直線方程也是在球坐標(biāo)系（r、θ、φ）下建立，并加入了射線仰角β的計(jì)算。令射線從T點(diǎn)傳播到A點(diǎn)，發(fā)射仰角為β0，方位角α（發(fā)射點(diǎn)正北向順時(shí)針偏移量），發(fā)射點(diǎn)所在地球的經(jīng)、緯度分別為φ0(0°≤φ0≤ 360°)、λ0(-9 0°≤ λ0≤ 90°，北緯為正，南緯為負(fù))，則射線的初始位置為：rT=OT， θT= π/2- λ0， φT=φ0，βT=β0。OT為地球半徑與初始點(diǎn)高度之和。射線以直線傳播，則到達(dá)A點(diǎn)（電離層底高）時(shí)的射線方程為：rA=OA，θA=θT- TOA ·co sα ， βA= π/2- TAO ， φA=φT+ TOA ·sinα / sinθT。OA為電離層底高，TOA和TAO可通過三角形正弦定理直接求得。

1.2 電離層段軌跡方程

在電離層AB段，以群路徑P'為自變量的射線方程可描述為[2-3]：

電離層段的射線方程涉及到群路徑步長(zhǎng)dP'的取值問題。步長(zhǎng)取的越小，結(jié)果越精確，但隨之相應(yīng)的問題是計(jì)算時(shí)間成倍的增加，因此采用變步長(zhǎng)技術(shù)來解決這一矛盾?？紤]到在射線追蹤中，影響射線傳播的主要因素是沿射線路徑的等離子體頻率的梯度，因此，根據(jù)等離子體頻率的梯度自適應(yīng)調(diào)整群路徑步長(zhǎng)。下面給出一種群路徑步長(zhǎng)的計(jì)算方式，如式（2）所示[4]：

其中，系數(shù)k的經(jīng)驗(yàn)值為0.01。只是給定一個(gè)k的取值還是不夠，因?yàn)楫?dāng)k取的較小時(shí)，雖然某些位置的步長(zhǎng)比較合適，但其他位置的步長(zhǎng)則可能會(huì)非常小，計(jì)算速度依然很慢；而當(dāng)k的取值較大時(shí)，其計(jì)算精度也無法保證。因此，可設(shè)定一個(gè)最小群時(shí)延步長(zhǎng)，使步長(zhǎng)不會(huì)小于該值，以保證計(jì)算速度；一個(gè)最大群時(shí)延步長(zhǎng)，使步長(zhǎng)不會(huì)大于該值，以保證計(jì)算精度[5]。

另外需要注意的是，射線在出電離層的時(shí)候要特殊處理，因?yàn)樯渚€在AB段的最后一步不可能正好落在B點(diǎn)，即r=OB（電離層底高）。此時(shí)應(yīng)根據(jù)dr = OB - rp （rp為上一步的高度）和式(3)重新確定dP'和其它位置參數(shù)的值，而B點(diǎn)的射線仰角 βB= arcsin(d r/ dP ')。

1.3 電離層底高到地面段軌跡方程

電離層底高到地面(BG)段的軌跡方程同樣根據(jù)直線幾何 關(guān) 系 直 接 計(jì) 算 ：rG=OG， θG=θB-BOG ·cosα ，φG=φB+ BOG ·si n α/ sinθB， βG= BOG + βB。OG為地球半徑，BOG依然可通過三角形正弦定理求得。需要說明的是，這里的βG為射線經(jīng)地面反射后的仰角，反射后的射線軌跡計(jì)算又與AB段的計(jì)算方法相同。

1.4 環(huán)球回波現(xiàn)象

射線經(jīng)電離層反射后，某些時(shí)候可能并不會(huì)到達(dá)地面，而是在地面與電離層之間來回反射，這種現(xiàn)象稱為環(huán)球回波現(xiàn)象[6]。如圖2所示。

作 B點(diǎn)至地球的切線 BG，由圖可知，當(dāng) OB· sin(βB+π/2)＞OG時(shí)，射線無法到達(dá)地面，發(fā)生環(huán)球回波現(xiàn)象。作O點(diǎn)至BC段的垂線OH，則OH即為BC段傳播過程中的最低高度，此時(shí)有：rH=OH，θH=θB-BOH ·cosα , βH=0,φH= φB+ BOH ·si n α/sin θB。OH、BOH可以由直線三角關(guān)系求得。至此，射線可能經(jīng)歷的各種情況都有了相應(yīng)的模式對(duì)其進(jìn)行計(jì)算。

2 計(jì)算結(jié)果比較

由以上算法可知，要實(shí)現(xiàn)電離層短波射線追蹤，必須知道電離層的等離子體頻率，等離子體頻率與電離層電子濃度和折射指數(shù)的關(guān)系為[7]：

式(3)是不考慮地磁場(chǎng)影響下電離層折射指數(shù)n的計(jì)算公式，Ne為電子濃度。因此只要知道電子濃度分布，就可以求解射線軌跡方程。

實(shí)際應(yīng)用中的射線軌跡顯示通常是建立在地面距離與高度的直角坐標(biāo)系上。因此這里就涉及到球坐標(biāo)系和直角坐標(biāo)系之間的轉(zhuǎn)換，地面距離x可通過下式計(jì)算：

設(shè)發(fā)射機(jī)位于地球表面的北緯35°、東經(jīng)113°位置，方位角0°，工作頻率15 MHz。電子濃度剖面利用亞大模型[8]計(jì)算。電離層參數(shù)：foE=4 MHz，h mE= 115 km， ymE= 20km，foF2 = 10 MHz， hmF2 = 330km， ymF2= 80km。這里給出發(fā)射仰角從6°～72°（遞增6°），k=0.01，= 0.1km，=1km條件下的快速算法和 dP'= 0.1km條件下的固定步長(zhǎng)算法的電離層短波射線軌跡比較圖。

從圖中可以看出，兩種算法的射線軌跡幾乎重合。同樣是利用MATLAB7.0 軟件進(jìn)行計(jì)算，使用快速算法用時(shí)53 s，而固定步長(zhǎng)算法則需245 s。

若改用附加赤道雙峰擾動(dòng)的準(zhǔn)拋物模型作為電離層模型[5]，設(shè)發(fā)射點(diǎn)位于地球表面的南緯 10°、東經(jīng) 113°位置，方位角20°，工作頻率10 MHz。電離層參數(shù)： fc0= 30 MHz，A1=0.3， A2=0.5，θ1=80°，θ2=100°， φ1=φ2= 10°，hb0=200，ym0=100，T= 10°， Crb = 0.9， C ym =-0 .4。射線仰角范圍與算法參數(shù)同前，得到等離子體頻率剖面以及快速算法（圖4中實(shí)線）和固定步長(zhǎng)算法（圖4中虛線）下的射線軌跡比較圖。

兩種算法的射線軌跡依然近乎重合，而快速算法用時(shí)28 s，固定步長(zhǎng)算法用時(shí)157 s。

這里之所以取固定步長(zhǎng) dP'= 0.1km進(jìn)行比較，是根據(jù)經(jīng)驗(yàn)認(rèn)為群路徑步長(zhǎng)取該值時(shí)即可滿足一般工程需要的準(zhǔn)確度[4]。

3 結(jié)語

快速算法在電離層以下采用球坐標(biāo)下的直線幾何近似，在電離層內(nèi)用球坐標(biāo)下的射線微分方程和變步長(zhǎng)技術(shù)，實(shí)現(xiàn)了短波射線軌跡的追蹤計(jì)算。該方法在保證計(jì)算精度的同時(shí)大大提高了計(jì)算速度，這對(duì)于需追蹤大量射線軌跡或?qū)?shí)時(shí)性要求較高的應(yīng)用系統(tǒng)來說更為方便實(shí)用。另外，由于所給出的射線方程將射線傳播過程中的高度變化、地心角變化以及方位角變化的因素都考慮在內(nèi)，因此快速算法也同樣適用于更為復(fù)雜的三維電離層電子濃度剖面。

[1] [1] 索玉成.電離層短波射線追蹤[J].空間科學(xué)學(xué)報(bào),1993, 13(04):306-312.

[2] [2] 郭杰,于大鵬,雷雪,等.基于數(shù)值射線追蹤的短波電離層傳播軌跡研究[J].通信技術(shù),2008,41(04):33-35.

[2] [3] 郭杰,于大鵬,雷雪,等.基于二維解析射線追蹤的短波電離層軌跡研究[J].通信技術(shù),2008,41(05):42-46.

[3] [4] 柳文,焦培南,王世凱,等.電離層短波三維射線追蹤及其應(yīng)用研究[J].電波科學(xué)學(xué)報(bào),2008,23(01):41-48.

[4] [5] 周向春,謝樹果,趙正予.變步長(zhǎng)技術(shù)在電離層射線追蹤中的應(yīng)用[J].武漢大學(xué)學(xué)報(bào),2001,47(03):355-358.

[5] [6] 宋錚,張建華,黃冶.天線與電波傳播[M].西安:西安電子科技大學(xué)出版社.

[6] [7] 種衍文,謝樹果,趙正予,等.一種新的射線追蹤方法-三區(qū)分劃處理法[J].電波科學(xué)學(xué)報(bào),2001,16(02):222-226.

[8] 江長(zhǎng)蔭,張明高,焦培南,等.-雷達(dá)電波傳播折射與衰減手冊(cè)[M].北京:國(guó)防科學(xué)技術(shù)工業(yè)委員會(huì),1997:30.