快速算法實現電離層短波射線追蹤

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

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

0 引言

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

1 快速算法

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

在TA和BG段，射線在電離層外，將射線看作直線傳播，利用幾何關系直接計算；而在AB段，則利用球坐標下的射線方程進行數值計算，并采用變步長技術，使算法達到準確且計算速度快的目的。下面將詳細分析各段的計算方法。

1.1 發射點到電離層底高段軌跡方程

由于電離層段使用的是球坐標下的射線方程，為了方便銜接計算，這里的直線方程也是在球坐標系（r、θ、φ）下建立，并加入了射線仰角β的計算。令射線從T點傳播到A點，發射仰角為β0，方位角α（發射點正北向順時針偏移量），發射點所在地球的經、緯度分別為φ0(0°≤φ0≤ 360°)、λ0(-9 0°≤ λ0≤ 90°，北緯為正，南緯為負)，則射線的初始位置為：rT=OT， θT= π/2- λ0， φT=φ0，βT=β0。OT為地球半徑與初始點高度之和。射線以直線傳播，則到達A點（電離層底高）時的射線方程為：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]：

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

其中，系數k的經驗值為0.01。只是給定一個k的取值還是不夠，因為當k取的較小時，雖然某些位置的步長比較合適，但其他位置的步長則可能會非常小，計算速度依然很慢；而當k的取值較大時，其計算精度也無法保證。因此，可設定一個最小群時延步長，使步長不會小于該值，以保證計算速度；一個最大群時延步長，使步長不會大于該值，以保證計算精度[5]。

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

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

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

1.4 環球回波現象

射線經電離層反射后，某些時候可能并不會到達地面，而是在地面與電離層之間來回反射，這種現象稱為環球回波現象[6]。如圖2所示。

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

2 計算結果比較

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

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

實際應用中的射線軌跡顯示通常是建立在地面距離與高度的直角坐標系上。因此這里就涉及到球坐標系和直角坐標系之間的轉換，地面距離x可通過下式計算：

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

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

若改用附加赤道雙峰擾動的準拋物模型作為電離層模型[5]，設發射點位于地球表面的南緯 10°、東經 113°位置，方位角20°，工作頻率10 MHz。電離層參數： 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。射線仰角范圍與算法參數同前，得到等離子體頻率剖面以及快速算法（圖4中實線）和固定步長算法（圖4中虛線）下的射線軌跡比較圖。

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

這里之所以取固定步長 dP'= 0.1km進行比較，是根據經驗認為群路徑步長取該值時即可滿足一般工程需要的準確度[4]。

3 結語

快速算法在電離層以下采用球坐標下的直線幾何近似，在電離層內用球坐標下的射線微分方程和變步長技術，實現了短波射線軌跡的追蹤計算。該方法在保證計算精度的同時大大提高了計算速度，這對于需追蹤大量射線軌跡或對實時性要求較高的應用系統來說更為方便實用。另外，由于所給出的射線方程將射線傳播過程中的高度變化、地心角變化以及方位角變化的因素都考慮在內，因此快速算法也同樣適用于更為復雜的三維電離層電子濃度剖面。

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

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

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

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

[4] [5] 周向春,謝樹果,趙正予.變步長技術在電離層射線追蹤中的應用[J].武漢大學學報,2001,47(03):355-358.

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

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

[8] 江長蔭,張明高,焦培南,等.-雷達電波傳播折射與衰減手冊[M].北京:國防科學技術工業委員會,1997:30.