一种基于地球同步轨道SAR的电离层电子密度反演方法

张 龙,田 野,李博骁

(1.中国电子科技集团公司，北京 100846；2.中国电子科学研究院，北京 100041)

0 引 言

电离层会对航天遥感、卫星通讯、导航定位等领域造成严重影响，因此利用现代技术手段准确获得实时电离层电子密度三维分布数据将具有非常重要的现实意义。同时，由于自然灾害发生前(如地震)以及发生时(如台风、火山爆发、海啸等)电离层电子密度将发生重大变化，因此准确获得实时电离层电子密度数据在地震预测、观测气象灾害发生时的电离层变化等方面具有重要的科研和应用价值。

电离层层析技术可准确反演出电离层电子密度空间分布，它基于观测区域内沿大量相互交叉的传播路径上的电离层总电子含量(Total Electronic Content: TEC)数值，通过层析反演技术获得一定空间范围内的电离层剖面结构。电离层层析技术由美国伊利诺伊大学Austen J于1986年首次提出[1]，Austen J指出结合卫星观测和层析计算可反演出二维空间的电离层电子密度分布。随后，在1988年Austen J对电离层层析理论进行了进一步详细阐述，并在1000 km卫星高度和地面分布3个和5个地面站假设的情况下实现了二维电离层层析反演[2]。

在电离层层析的电子密度计算算法中，提出最早的为代数重构算法(ART)[3]。该算法是在电子密度重建过程中逐步改善电子密度迭代初始值，使结果逐渐趋近于最佳电子密度结果的过程。ART算法的最大问题在于在某些条件下迭代出的电子密度数值会为负值，不符合实际情况，因此一种改进的迭代算法应运而生。乘法代数重构算法(MART)[3]是目前最常用的一种电离层电子密度迭代反演算法。MART迭代来自于熵最优算法，在利用该算法进行电离层三维层析反演时，先对各个像素网格设定一个电子密度初值，然后根据初始值通过迭代求出一次近似值，再通过迭代求出二次近似值，如此继续，直到迭代满足需求阈值为止。与ART算法相比，该算法可保证迭代结果为正值，同时具有收敛速度快等优点，因此它已成为电离层层析算法中最常用的一种迭代算法[4-7]。其它电子密度计算算法还有Fremouw提出的离散反演(DIT)算法[8]，Na H和Lee H提出的正交分解算法[9]等。

目前的主流电离层层析方式是基于GPS卫星，它利用布设在地面的观测站或电离层监测仪接收GPS卫星信号，获得GPS信号传播路径上的电离层TEC数据，并基于获得的电离层TEC数据进行电离层层析计算。由于GPS卫星可实现全球覆盖，因此能够保证在地球上任意位置任意时刻观测到至少四颗GPS卫星；此外，国际大地测量协会(IAG)在全球布设了几百个电离层观测站，且数目仍在不断增加，这为电离层TEC数据获取提供了可靠的保证。因此，目前基于GPS卫星的电离层层析技术取得了长足发展和广泛应用。

然而，基于GPS卫星的电离层层析方法也具有非常明显的缺点。由于电离层观测站无法精确覆盖全球任意角落，因此无法提供任意位置的精确电离层TEC数据。此外，由于部分地区自然环境条件的限制，无法在恶劣环境下布设电离层监测仪来获取GPS路径上的TEC数据。因此，基于GPS卫星的电离层层析方法具有一定局限性。

近几年出现的基于星载SAR卫星的电离层层析方法可有效解决该问题。星载SAR信号在传播过程中会受到电离层效应的影响，因此SAR信号中会携带明显的电离层信息。通过对SAR信号进行处理可反演出信号传播路径上的电离层TEC数据，同时，通过对成像区域内各条射线路径上的SAR信号进行处理可反演出三维空间中的电子密度数据。与传统电离层TEC测量和电子密度反演方法比，基于星载SAR的电离层参数反演方法可实现电离层参数的实时估计和反演，这对星载SAR信号处理具有重要意义。但是，基于星载SAR的电离层层析方法均是基于低轨星载SAR数据，它通过接收地面强散射点的星载SAR回波信号，从信号中反演出星载SAR信号路径上的电离层TEC数据，并基于该数据进行电离层层析反演[10,11]。基于星载SAR的电离层层析几何结构和层析结果示意图如图1所示。

目前，基于星载SAR的电离层层析方法均是利用低轨星载SAR数据，由于低轨星载SAR的重访时间长、覆盖范围小，因此无法对特定区域进行实时大范围的电子密度反演。因此，亟需研究新的电离层电子密度层析方法，以实现对特定区域的三维空间实时全覆盖层析反演。

图1 基于星载SAR的电离层层析几何结构和层析反演结果[11]

此外，在电离层TEC反演过程中，既有的基于星载SAR数据的电离层TEC反演方法较少考虑电离层延迟误差以外其它误差的影响。地球同步轨道SAR(GEO SAR)信号能够完整地穿越大气层，并且由于其长孔径时间和大覆盖范围的特性导致波束覆盖区域内空间环境复杂，因此需考虑电离层延迟误差外其余误差的影响，如对流层误差、电离层多重散射误差、随机噪声误差等[12]。

本文提出了一种基于地球同步轨道SAR数据的电离层三维电子密度反演方法，可实现对特定大范围区域的实时连续高精度三维电子密度反演。文章第二节首先提出了一种复杂环境下的电离层TEC参数反演方法，综合考虑背景电离层、电离层多重散射、对流层等环境误差影响下的电离层TEC反演方法。第三节利用地球同步轨道SAR回波数据和星地三维运动关系，提出了一种基于GEO SAR的电离层电子密度反演方法。在第四节，利用电离层实测数据进行了TEC数据反演和三维空间电子密度反演，结果表明了使用方法的有效性。

1 电离层TEC参数反演

根据既有研究，GEO SAR信号在穿越大气层的过程中会受到电离层、对流层等多重因素影响，由此引入的相位误差为[12]

(1)

在式(1)中，第一项为电离层延迟误差，第二项为电离层多重散射误差，第三项为对流层误差，第四项为随机噪声误差。f0为信号载频，fr为信号距离向频率，Ciono为电离层多重散射系数，ΔRtroposphere为对流层引起的信号延迟，Δφrandom为随机相位误差。

在传统低轨SAR中，由于合成孔径时间非常短，因此在孔径时间内式(1)中的气象参数不随时间变化，即在孔径时间内可将式(1)中的气象参数视为定值。然而，由于GEO SAR存在多达几百秒至上千秒的孔径时间，因此需要考虑在孔径时间内气象参数随时间的变化情况对GEO SAR成像的影响。考虑气象参数的变化，式(1)可表示为

(2)

在式(2)中，随时间变化的气象参数会导致出现随时间变化的相位误差，该误差将对GEO SAR成像造成影响。根据本文前面章节的分析，可用高阶多项式描述电离层TEC和对流层延迟误差，TEC和对流层延迟误差可分别表示成下式

(3)

ΔRtroposphere(ta)=ΔRtroposphere_0+q1·ta+

(4)

在式(3)和(4)中，TEC0和ΔRtroposphere_0分别为TEC(ta)和ΔRtroposphere(ta)的常数部分，ki,i=1…n和qi,i=1…n分别为TEC(ta)和ΔRtroposphere(ta)对时间的各阶导数。

由于电离层多重散射系数的成因和变化形式较为复杂，且尚未有实时电离层多重散射系数的实测数据，因此尚不清楚它随时间的变化形式能否用高阶多项式形式描述。因此，在这里我们先认为它在合成孔径时间内为一定值[11]，关于电离层多重散射随时间变化的数学建模和反演方法将在未来进行具体研究。

至此，我们已经建立了大气层环境影响下的GEO SAR信号模型，如式(2)所示。接下来将研究GEO SAR电离层TEC反演方法，即研究如何从受大气环境影响的GEO SAR回波数据中反演出电离层TEC数据。

基于GEO SAR回波数据反演电离层TEC具体处理流程为：首先，接收地面PS点处的GEO SAR回波信号；然后，为了获取随合成孔径时间变化的TEC数据，需要对GEO SAR数据进行子孔径划分，并假设在每一个子孔径时间内TEC为一恒定值；其次，为了获取每个子孔径内的电离层TEC数据，对子孔径内的GEO SAR回波信号进行分谱处理，将其划分成若干个子带信号；紧接着，对各子带GEO SAR信号分别进行保相成像处理，获取各个子带信号的残余相位误差，并解算出子孔径内的电离层TEC数值；最后，对各个子孔径内获得的TEC数据进行方位向最小二乘拟合，获得随孔径时间变化的电离层TEC数据。整体流程图如图2所示。

图2 基于GEO SAR回波数据反演电离层TEC流程图

根据电离层TEC反演流程可以看出，在电离层TEC的反演过程中有三个关键步骤：一是方位向划分子孔径，二是距离向分谱，三是根据残余相位误差解算电离层TEC。下面分别给出这三个关键步骤的具体处理方法。

1.1 划分子孔径

为了获得随时间变化的电离层TEC数据，我们需要将GEO SAR的超长合成孔径时间划分成若干个子孔径时间，在每个子孔径内，将电离层TEC视为一恒定值。为了确保电离层TEC反演精度，因此确定子孔径时间长度的原则是在每个子孔径时间内电离层TEC的变化对聚焦的影响可以忽略。

引起GEO SAR图像散焦的主要是TEC随时间的二阶变化率，电离层引起的方位向二次相位误差为[12-13]

(5)

根据SAR成像原理，当式(5)所示的二次相位误差小于π/4时，电离层不会对GEO SAR聚焦造成影响。因此，可得最大子孔径时间为

(6)

1.2 分谱处理

在方位向子孔径划分完成后，需要对距离向进行分谱处理，这主要是基于两个原因。首先，从式(2)可以看出，大气相位误差表达式中存在三个未知参量：电离层TEC、多重散射系数、对流层延迟，此外还有随机噪声的影响。根据图2所示的电离层TEC反演流程，在对GEO SAR回波信号进行保相成像处理后，可提取出成像后的残余相位误差，该误差即为式(2)所示的大气相位误差。然而，由于该相位误差中存在三个未知参数，因此无法通过一个方程求解出三个未知参数。对距离向进行分谱处理后，可以获得不同信号频率对应的相位误差表达式，为后续的电离层TEC解算做铺垫。其次，由于成像后的残余相位误差存在相位模糊现象，即读取到的残余相位误差和真实相位误差之间相差2π的整数倍，无法反应真实的相位信息，也无法应用它解算电离层TEC。在进行距离向分谱之后，相邻子带残余相位之差小于2π，此时将不存在相位模糊现象。联立相邻子带残余相位之差，即可建立方程组解算出电离层TEC数据。

将带宽为B的回波信号划分成N个子带，每个子带的带宽为Bsub=B/N。子带划分示意图如图3所示。假设全孔径GEO SAR回波信号被划分为M个子孔径，第m个子孔径时刻第n个子带处的信号相位误差Δφ(m,n)可以表示为

(7)

其中，f0n为第n个子带的信号中心频率，Δφrandom_n为第n个子带的随机噪声误差，TECm、Ciono_m和ΔROthers_m分别为第m个子孔径处的TEC、电离层多重散射系数和对流层延迟数据。

图3 分谱示意图

为了消除相位模糊现象的影响，我们对各个相邻子带的残余相位做差，可得

(8)

为方便后续分析，将未知气象参数均用符号进行表示，上式可简化为如下形式

Δφm,n=an·TECm+bn·Ciono_m+cn·ΔROthers_m+dn

(9)

其中

Δφm,n=Δφ(m,n+1)-Δφ(m,n)

(10)

(11)

(12)

(13)

dn=Δφrandom_n+1-Δφrandom_n

(14)

在后续分析中，我们将利用式(9)～(14)解算电离层TEC数据。

然而，在实际操作中，子带带宽过宽或过窄都会影响最终TEC反演质量。限制子带带宽的主要有两个因素：一是，子带带宽上限受到相位模糊现象约束，要求相邻子带相位差不会再次出现相位模糊现象，即相邻子带相位差的绝对值不应超过2π；二是，子带带宽下限受到子带回波信号信噪比限制，子带带宽过小会导致子带回波信号信噪比过低，影响反演精度。下面分别就这两点进行具体分析。

子带带宽上限应保证相邻子带相位差的绝对值不超过2π，对于L波段GEO SAR系统，不同TEC大小的情况下，相邻子带的相位差如图4所示。

图4 子带带宽、TEC和相邻子带相位差的对应关系

考虑到在正常情况下，电离层TEC一般不会超过60TECU。通过图4的仿真分析可以看出，为了不失一般性，我们通常选择子带带宽不超过10 MHz。

根据雷达方程，信噪比表达式为

(15)

其中，Pt为平均发射功率;G为天线增益;A为天线面积;ηT为损耗因子;σ0为后向散射系数;θ为入射角;ρr为距离向分辨率;N0为噪声功率谱;vs为平均速度;R为GEO SAR卫星到目标的速度。

从式(15)可以看出，子带数量越多，每个子带信号对应的成像结果距离向分辨率越小，相应的信噪比也越低，会降低相位误差提取精度，并进而会降低TEC反演精度。设GEO SAR信号带宽为40 MHz，图像中PS点信噪比为30 dB，在距离向分别划分子带数量为4个、5个、6个、7个、8个、9个，对应的各子带信号PS点成像结果的信噪比分别为18 dB、16 dB、14.5 dB、13 dB、12 dB、11 dB。相应的以均方根误差衡量的TEC反演精度为0.3TECU、0.9TECU、1.5TECU、4.5TECU、13.8TECU、22.1TECU。

从图5的仿真结果可以看出，信噪比越差，TEC反演精度也越差。因此，在实际操作时，应尽量少划分子带，同时要选择信噪比尽量大的地面PS点目标。

图5 信噪比，子带数量与TEC反演精度的关系

1.3 解算电离层TEC

在进行完方位向划分子孔径和距离向划分子带之后，提取各子带对应的残余相位误差并将其做差，可得第m个子孔径处的相位误差方程为

(16)

其中，ai,bi,ci,di,Δφm,i分别如式(10)～(14)所示。为了分析表述方便，上述方程组可以记作

A·x+ε=φ

(17)

其中

(18)

式(17)所示的方程组为超定方程组，并且存在由随机噪声引入的相位误差ε。由于噪声是随机的，无法获取其具体数值或表达式，在后续公式推导时将其忽略。因此，可将式(17)视为

A·x=φ

(19)

下面将针对式(19)所示的超定方程进行求解。求解的第一步首先需要将其化为正定方程组，方法为

ATA·x=ATφ

(20)

因此，式(20)可化为

M·x=B

(21)

由此可得式(16)中各未知参数的通解为

(22)

(23)

(24)

式(22)～(24)即为各个参数的最优解。

至此，我们已经获得了在第m个子孔径时间内的TEC、电离层多重散射系数和对流层延迟的最优解，分别对所有子孔径内的GEO SAR回波数据使用上述方法求解出各个子孔径内的TEC、电离层多重散射系数和对流层延迟数据，并将各个子孔径内的数据沿合成孔径时间进行高阶多项式拟合，即可获得随时间变化的电离层TEC和对流层延迟数据。

2 基于地球同步轨道SAR的电离层三维电子密度反演

与传统GPS电离层层析和低轨星载SAR电离层层析相比，GEO SAR层析具有时间和空间覆盖范围大、操作简便、实时性强等突出优点。

图6所示为单颗GEO SAR卫星和ERS-1卫星在24小时时间内进行三维层析时射线在反演空间区域中的空间占比(即有射线穿过的空间范围占总空间范围的比重)。从图6可以看出，GEO SAR卫星波束在大部分时间内能够覆盖超过50%的反演区域，而利用ERS-1卫星进行电离层三维层析反演时只有极少部分时间可以照射到少量反演区域，而在超过95%的时间里都无法照射到反演区域。由此可见GEO SAR层析具有非常大的优势。

图6 单颗GEO SAR卫星与ERS-1卫星电离层三维层析时空占比示意图

虽然与低轨星载SAR层析相比，GEO SAR层析的优势已经非常明显，然而使用单颗GEO SAR卫星进行电离层层析反演尚存在无法对层析反演区域进行连续观测的缺点。因此在后续研究中可以通过进行合理的GEO SAR卫星编队设计以弥补单颗GEO SAR卫星层析时的缺点。

GEO SAR电离层层析几何结构如图7所示，O-XYZ为地固坐标系，GEO SAR卫星坐标为S=(Xs(n),Ys(n),Zs(n))，其中Xs(n),Ys(n),Zs(n)分别表示GEO SAR卫星在第n时刻的三维坐标值。地面分布有若干个PS点，PS点坐标为Pi=(Xip,Yip,Zip)，Pi表示第i个PS点的坐标，Xip,Yip,Zip分别表示第i个WTPS点的三维坐标值。在进行电离层层析反演时，需要将三维空间划分为若干个子网格，在每个子网格内认为电离层电子密度为一定值。

图7 GEO SAR电离层三维层析几何结构示意图

考虑单颗GEO SAR卫星依然有接近一半的时间无法照射到反演区域，这里我们设置4颗GEO SAR卫星组成卫星编队进行电离层层析。GEO SAR卫星编队如图8(a)所示，4颗卫星星下点轨迹等间隔分布。此时，电离层层析的时空分辨率如图 8(b)所示，可见4颗GEO SAR卫星组成的编队可以实现任意时间对接近100%的层析空间的全覆盖。同时，由于GEO SAR卫星可以从不同方向照射层析区域，因此在层析计算时可以获得更多的电离层信息，进而能够提高电离层层析精度。

基于GEO SAR的电离层层析流程如图9所示。当GEO SAR卫星向地面发射电磁波时，通过接收地面PS点处的回波信号并利用上文给出的电离层TEC反演方法，从GEO SAR回波信号中提取出信号路径上的TEC数据用于电离层层析反演。同时，将待反演区域划分为一系列空间网格，根据各个时刻GEO SAR卫星和PS点的位置关系以及网格划分情况，计算出各个空间网格中各条射线穿过的截距长度，建立投影矩阵。最后，利用国际参考电离层(IRI)[13-14]提供的三维空间中电离层电子密度数值作为电离层层析反演的初始值数据，并根据MART迭代算法计算出三维空间中的精确电子密度分布情况。

图8 GEO SAR四颗卫星编队及电离层层析时空覆盖率

图9 使用GEO SAR回波数据进行电离层层析反演流程图

根据图7的几何关系可知，在第n时刻从GEO SAR卫星到第i个PS点的射线方程可表示为

(25)

其中，X,Y,Z表示该射线与层析网格面的交点坐标;K表示射线与网格面的交点在射线上的位置。

在进行层析网格划分时，通常以地理经度、纬度和海拔高度为参考进行划分，分别对应层析网格面中的经度面、纬度面和高度面。定义高度面是以地心为球心的球面，半径为反演区域的平均地球半径R加上该高度面到地面的高度H

X2+Y2+Z2=(R+Hm)2

(26)

其中，Hm表示第m个高度面对应的地面高度。

定义经度面为一个穿过Z轴的平面，且与起始子午面夹角为经度J，其方程为

tanJp·X-Y=0

(27)

其中，Jp为第p个经度面对应的经度值。

定义纬度面为一个由某点与地心连线绕Z轴旋转所形成的圆锥面，该面上任一点的纬度都是相等的，若纬度为W，其方程为

Z2=tan2Wq·(X2+Y2)

(28)

其中，Wq表示第q个纬度面对应的纬度值。

将式(25)分别与式(26)(27)(28)联立，可求解出射线与各个高度面、经度面、纬度面的交点，即可计算出射线在各个三维层析网格中的截距长度。同时，使用MART迭代算法即可求解出三维空间中精确的电子密度数据。

在第k次迭代中，MART算法的修正公式为

(29)

3 实测数据仿真

这部分我们将对GEO SAR电离层三维层析进行仿真分析。根据前一节的分析，为了实现实时全空间的电离层层析覆盖，我们设置4颗GEO SAR卫星组成卫星编队，同时假设地面分布若干PS点。根据各PS点处的GEO SAR回波信号，反演出各GEO SAR卫星到各PS点路径上的TEC数据。GEO SAR卫星编队进行电离层层析示意图如图10所示。

图10 GEO SAR卫星编队电离层层析示意图

仿真分析时，利用2015年5月13日UT8:00的国际参考电离层(IRI)数据[14]作为三维空间中电子密度分布的初始数据，在进行电离层层析计算时将其作为初始迭代值。同时，利用2015年5月17日UT8:00的IRI数据作为三维空间中电子密度分布的目标精确值，将电离层层析的结果与该目标值进行比较，以验证电离层层析的精度。初始迭代值和目标精确值之间存在大约20%的误差。电子密度初始迭代值和精确值均位于100°E～120°E、10°N～30°N、高度200～600 km的空间范围。在层析计算时，将该空间范围进行网格划分，空间网格间隔在纬度维为0.5°，经度维为1°，高度维为40 km。

图11显示了电子密度初始值、迭代值和精确值之差，从图11 (a)可以看出，在没有进行电离层层析计算时，电子密度初始值和电子密度精确值之间具有较大误差。从图11(b)可以看出，在进行完电离层层析运算之后，电子密度迭代值和目标精确值之间的误差已经大大缩小，此时迭代值已经非常接近目标精确值。

为了更好的显示电离层层析效果，分别取出110°E、17°N和110°E、20°N处的电离层层析结果。从图12(a)(b)的电离层层析反演结果可以看出，使用GEO SAR卫星可以实现良好的层析效果，初始值与精确值之间存在约20%的误差，但经过电离层层析之后的反演值已经非常逼近层析精确值。

图11 电子密度初始值、迭代值和目标精确值之间误差(110°E，10°～30°N)

图12 电子密度初值、迭代值和精确值电子密度剖面图

另外，由于GEO SAR电离层层析时使用了卫星编队，因此它具有多角度层析的特点，即比单颗GEO SAR卫星具有更加丰富的电离层射线信息，从而会有更加精确的电离层层析结果。为了比较GEO SAR多角度层析对反演精度的影响，我们仿真了单颗GEO SAR卫星的层析结果，并将其与4颗GEOSAR卫星组成的编队层析结果进行比较，比较结果如表1所示。从中可以看出，使用GEO SAR卫星编队进行电离层层析可使得层析精度提高约20%。

表1 单颗GEO SAR卫星与4颗GEO SAR卫星层析精度比较

4 结 语

文章提出了一种基于GEO SAR回波数据的电离层电子密度高精度反演方法。首先，在穿过大气层的GEO SAR回波信号模型的基础上，提出了一种基于GEO SAR回波信号的实时时空变电离层TEC数据反演方法，它利用距离向划分子带、方位向划分子孔径和最小二乘原理，基于GEO SAR回波信号实现了实时电离层TEC高精度反演。然后，利用GEO SAR卫星并结合电离层层析理论，研究了GEO SAR电离层三维层析的概念、几何构型、卫星编队以及实现流程，提出了一种基于GEO SAR数据的全新的实时三维空间电子密度反演方法。最后，基于IRI实测数据实现了三维空间中实时高精度的电离层电子密度反演。