快速算法实现电离层短波射线追踪

孙 方， 康士峰， 赵振维， 王红光

(中国电波传播研究所青岛分所，山东 青岛 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.