基于星载AIS 数据的TEC 测量方法

2022-06-20罗坳柏龙燕飞陈利虎方涵先余孙全倪久顺

空间科学学报 2022年3期

罗坳柏龙燕飞陈利虎方涵先余孙全倪久顺

1(国防科技大学空天科学学院长沙 410073)

2(国防科技大学气象海洋学院长沙 410073)

0 引言

基于星载船舶的自动识别系统（AIS）采用开放式广播技术，周期性播发自身位置、速度等信息，广泛分布于全球。近年来，星载AIS 发展迅速，目前能够利用卫星平台在全球范围内实时接收信号。AIS 信号从船舶传播到卫星时会穿过地球大气层。0.1～12 GHz 频率范围的电磁波通过电离层时，可能引起信号特性改变。线极化波在电离层中传播时，由于等离子介质存在各向异性，受到地磁场的作用，极化平面将发生旋转，这种现象称为法拉第效应。

电离层电子总含量TEC 是描述电离层形态和结构的主要参量之一，对通信、导航、地震预测等电离层相关技术具有重要影响。目前常用的TEC 测量平台有星载广播式自动相关监视技术ADS-B（Automatic Dependent Surveillance-Broadcast）、星载合成孔径雷达SAR（Synthetic Aperture Radar）和全球定位系统GPS（Global Positioning System）。常用的TEC测量理论有法拉第旋转法、差分群时延法和差分多普勒法。文献[1]以GPS 卫星为测量平台，提出了一种双频时延与相位联合的TEC 测量方法。在此基础上，文献[2]进一步提出三频时延与相位联合的TEC 测量方法，提高了相位积分常数的测量精度。文献[3,4]以6 颗卫星组成的COSMIC 星座为测量平台，采用GPS 掩星技术进行TEC 测量。但是基于GPS 建立的电离层模型的时间和空间分辨率都相对较差，时间分辨率一般为2 h，空间分辨率在经纬度方向上分别为5°和2.5°[5]。近年来，低轨卫星在提高电离层模型精度方面发挥了重要作用[6]。文献[7]使用海军卫星导航系统NNSS（Navy Navigation Satellite System）数据进行TEC 测量。文献[8–10]证明了利用低频SAR 数据进行TEC 测量的可行性。文献[11,12]以星载ADS-B 为测量平台，采用射线追踪技术确定单个电磁波的特性，然后计算法拉第旋转角，并将其转换为沿射线路径的TEC。文献[13]证明了AIS 信号可以作为TEC 测量的数据源。

以上提到的TEC 测量方法需要在地面布置相关基站，这将消耗大量人力物力，并且基站存在难以在海洋、极地和荒漠等地区部署的难题。本文研究认为当使用两个相互垂直的AIS 线极化天线接收同一个AIS 信号时，接收到的信号功率是不同的，接收功率大小与电磁波极化面的位置有关。具体研究过程：通过两台AIS 接收机接收功率之间的关系，计算出AIS 信号发生法拉第旋转之后的角度，根据船舶和卫星上天线的安装位置、卫星的速度和位置信息、船舶的位置信息，计算AIS 信号发生法拉第旋转之前的角度，通过这两个角度计算出法拉第旋转角，进而根据法拉第旋转角计算TEC。此外，分析了接收机硬件误差和观测参数误差对结果造成的影响，并选取天拓五号卫星上的AIS 数据进行实验验证。

1 方法原理

用于求解法拉第旋转角的卫星AIS 载荷主要包括两台AIS 接收机和两个垂直安装的AIS 线极化天线，如图1 所示。其中AIS 一体化接收机由高灵敏度AIS 接收模块和电源组件等组成。AIS 天线用于接收频率为156～163 MHz 的AIS 信号。AIS 接收机可以完成对AIS 信号的采样、下变频、AD 采样、解调和组包，再统一传送给数传分系统，当卫星过境时，将AIS 报文（包含时间戳、报文功率、目标原始位置）传输给地面测控站，最后由地面测控站将解析后的AIS 数据传输给用户。

如图2 所示，参考天拓五号卫星上AIS 天线的安装位置，以卫星质心为原点，建立直角坐标系，记为坐标系1。其中+z1指向地心，+x1指向卫星飞行方向，y1轴的方向遵循右手定则。两个AIS 天线分别安装在±x1面上，与z1轴成45°安装。其中AIS 天线1 在–x1面，AIS 天线2 在+x1面上。

使用卫星工具包STK（Satellite Tool Kit）得到的卫星位置和速度以及从AIS 报文中读取的船舶相关位置坐标是在WGS-84 坐标系下显示的。由于下文的计算基本都是在地心地固坐标系下进行的，因此这里需将这些信息转换在图2 中的地心地固坐标系下表示。以WGS-84 椭球为基准，空间中某一点的坐标为(λ,γ,H)，通过下式可以变换到地心地固坐标系(x2,y2,z2)中：

图1 用于TEC 测量的星载AIS 系统Fig.1 Schematic diagram of the space-based AIS system for TEC measurement

图2 用于TEC 测量的相关坐标系Fig.2 Relevant coordinate system used for TEC measurement

式中，N表示椭圆球卯酉圈的曲率半径，e表示偏心率，a表示点所在轨道半长轴，b表示点所在轨道半短轴。

如图3 所示，为便于进行计算过程描述，将发生法拉第旋转之前的极化面记为平面1，发生法拉第旋转之后的极化面记为平面2，两个AIS 天线所在的平面记为平面3，过AIS 天线2 的方向向量并且与平面3 垂直的平面记为平面4。该方法的工作流程如图4所示。

图3 基于星载AIS 数据的TEC 测量方法Fig.3 TEC measurement method based on space-based AIS data

图4 基于星载AIS 数据的TEC 测量方法流程Fig.4 TEC measurement method based on spacebased AIS data

1.1 AIS 报文选取

选取图1 中接收机A 和B 收到的同一条AIS 报文，可以得到船舶的经纬度信息。船舶上AIS 天线安装位置始终垂直于船身表面。卫星的位置和速度信息可通过查询卫星星历得到。

假设在地心地固坐标系下，卫星的位置坐标为(x,y,z)，卫星速度的方向向量坐标为v=(x˙,y˙,z˙)，船舶的位置坐标为(x0,y0,z0)。根据这些条件可以进一步得出：船舶上对天的AIS 天线的方向向量m=(x0,y0,z0)，卫星到地心的方向向量s=(−x,−y,−z)，船舶到卫星的方向向量f=(x −x0,y −y0,z −z0)。由空间几何关系可知，平面1 的法向量为

平面3 的法向量为

根据式（6）和式（7），可以得到平面1 与平面3 的夹角

理想条件下，AIS 信号沿船舶到卫星的方向直线传播，可以得到平面3 的法向量n与AIS 信号传播方向的向量f的夹角

理想条件下，船舶与卫星的位置关系如图3 所示。根据式（8），判断平面1 与平面3 是否近似垂直，这里规定二者夹角一般与90°相差不超过10°，记为条件1。根据式（9），判断AIS 信号的传播方向与平面3 的法向量是否近似平行。由于卫星与船舶存在一定的高度差，因此这里规定二者夹角一般与0°相差不超过30°，记为条件2。如果以上两个条件都满足，利用选中的AIS 报文中的功率和位置信息进行法拉第旋转角计算；如果不满足，则重新寻找同时满足以上两个条件的报文。

1.2 法拉第旋转角计算

根据1.1 节选中的这组两台接收机同时收到的同一条AIS 报文，记录其功率信息，可以计算得出平面2 与平面4 的夹角

式中，PA为AIS 接收机A 的接收功率值（单位dBm），DA为天线1 的接收增益，∆PA为接收机A 存在的功率误差，PB为AIS 接收机B 的接收功率值（单位dBm），DB为天线2 的接收增益，∆PB为接收机B 上存在的功率误差。

因为发射源为同一个，发射功率、发射天线增益和传输损耗都是一样的，所以式（10）中的天线增益只需减去接收天线的增益。船舶位置可通过查询AIS报文得知，卫星位置通过查询卫星星历得知，卫星上AIS 天线的方向图和安装位置事先已知，根据这些信息可以计算得出天线方向图中的方位角和仰角，进而查收接收天线的增益。

如图2 所示，为便于查找AIS 天线的方向图，这里的仰角指入射波与坐标系1 中+z1轴之间的夹角。+z1轴的方向向量为卫星到地心的方向向量s=(−x,−y,−z)，入射波的方向为船舶到卫星的方向向量f=(x −x0,y −y0,z −z0)，可得到入射波相对于接收天线的仰角

如图2 所示，为便于查找AIS 天线的方向图，这里的方位角指的是坐标系1 中入射波在Ox1y1平面上的投影向量相对于+x1轴偏转的角度。+x1轴指向的方向为卫星速度方向v=设s′为卫星到地心的方向向量s的单位向量，f′为船舶到卫星的方向向量f的单位向量。入射波的方向向量为f′，该向量在Ox1y1平面上的投影向量为f0=cosφs′+f′。入射波相对于接收天线的方位角

设v′为卫星速度的方向向量v的单位向量。AIS天线1 的方向向量

AIS 天线1 的方向向量即为平面4 的法向量。根据式（6）和式（13），可以得到平面1 与平面4 的夹角

根据式（10）和式（14）可以得到AIS 信号穿过电离层时发生的法拉第旋转角

1.3 TEC 计算

如图5 所示，假设电离层为单层薄壳模型，可以把接收机和卫星连线与电离层薄层的交点称为电离层穿刺点IPP（Ionospheric Pierce Point）[14]。图5 中R为地球的平均半径，一般取6371 km，天拓五号卫星的高度约为500 km；H为电离层薄壳的高度，一般为300～450 km[15]。

图5 基于电离层薄层模型计算穿刺点Fig.5 Calculating the puncture point based on the ionospheric single layer model

(θel,θAz)代表卫星相对于船舶的高度角和方位角坐标。根据船舶的位置坐标(x0,y0,z0)和船舶到卫星的方向向量f=(x −x0,y −y0,z −z0)，可以求出高度角为

式中，船舶的位置坐标用船舶上对天的AIS 天线的方向向量m表示。(φlat,φlon)为船舶的纬度和经度坐标，(ψlat,ψlon)为卫星的纬度和经度坐标。方位角为

z′为IPP 处的天顶角，有

为求出IPP 处垂直方向的电离层TEC，需先知道IPP 处的位置信息。IPP 处纬度和经度坐标(φlat,φlon)可由下式求出：

式中，

TEC 可以写为[16]

法拉第旋转角与TEC 的关系如下：

式中，l表示射线路径，ne表示传播路径l上每一点的电子密度，bz表示传播路径l上每一点磁场强度在该方向上的分量，e=1.6021892×10–19C 表示电子的电荷，c=3×105km·s–1表示真空中的电磁波传播速度，me=9.10956×10–31kg 表示电子的质量，ε0=8.854187817×10–12C2·N–1表示真空中的介电常数，w为AIS 信号的角频率。

联立式（22）与式（21），并代入已知量可得

由于采用电离层单层模型，式（22）中的bz近似为IPP 处的磁场强度Bavg。这里Bavg指平均电磁场强度。目前已知IPP 处的纬度和经度坐标为 (φlat,φlon)，Bavg可根据国际地磁场参考模型IGRF（International Geomagnetic Reference Field）[17]计算得出。

2 误差分析

在使用本方法进行TEC 测量的过程中，存在许多对测量结果有影响的误差来源，误差分为硬件设备误差和观测参数误差。对这些误差影响进行分析有助于提高测量精度。目前发现的硬件设备误差主要有AIS 接收机报文功率的读数误差及测量误差，观测参数误差主要有卫星姿态变化误差、船舶姿态变化误差、船舶位置误差、卫星星历误差以及磁场强度Bavg误差。

AIS 接收机的报文功率读数误差。AIS 报文中的16 进制功率值对应的实际功率值精度为1 dBm，造成的功率误差最大可达0.5 dBm。

AIS 接收机的报文功率测量误差。实际利用AIS 接收机对模拟器报文进行功率标定，信号输入为恒定时，功率差变化范围约±10%。接收机的报文功率测量误差包括两接收机射频通道的放大倍数误差（固定值可标定去除）、AD 量化位数误差、数字滤波器、下变频和FFT 等误差。由功率差造成的角度θ1的测量误差与角度值相关。例如在90°和0°附近，功率的微小偏差会造成角度的剧烈变化。以约45°为例，功率测量值浮动10%，幅度变化3.3%，角度测量误差范围为±6.6%。

卫星姿态的变化误差。卫星姿态趋于稳定，姿态角误差约0.1°左右，造成的法拉第旋转角误差可以忽略不计。

船舶姿态的变化误差。选取的船舶一般处于平稳航行状态，航行过程中船身摆动一般也很小，约为1°。船身摆动可能造成天线摆动，使得极化面发生偏转，造成的法拉第旋转角最大误差约为1°。

船舶的位置误差。AIS 报文自带GPS 坐标，其中单条报文坐标精度优于100 m，如果奇/偶报文联合解析，坐标精度优于20 m。由于船舶与卫星的相对距离为100 km 量级，因此造成的法拉第旋转角误差可以忽略不计。

卫星星历误差。卫星星历误差中卫星的位置和速度变化误差为100 m 量级。卫星的速度约为7.5 km·s–1，船舶与卫星的相对距离为100 km 量级，卫星位置及速度的变化误差造成的法拉第旋转角误差可以忽略不计。

磁场强度Bavg误差。虽然电离层有几百千米厚度，但是电离层的大多数电子实际集中在一个相对薄的层中，因此式（22）中的地磁场强度分量bz可以简化为薄层中的平均地磁场强度分量Bavg。磁场强度Bavg误差大小与电离层薄层高度的选取有关。一般情况下，电离层薄层高度的变化会引起最大3.2°的穿刺点经纬度差异和最大9.62%的电离层模型精度误差，进而引起磁场强度Bavg变化。

3 实验验证

选取天拓五号卫星上的AIS 数据进行实验验证。2020 年8 月23 日卫星上接收到的AIS 数据如图6 所示。

图6 根据2020 年8 月23 日天拓五号卫星接收到的AIS 数据得到的船舶分布位置Fig.6 Distribution of ships based on the AIS data received on the Tiantuo V satellite on 23 August 2020

利用STK 载入天拓五号卫星的星历获取2020年8 月23 日05:00 UT 卫星的位置和速度信息。卫星经度为121.157°W，纬度为15.495°S，高度为497.811 km，纬度变化率为0.063229 (°)·s–1，经度变化率为0.013092 (°)·s–1，高度变化率为–0.006980 km·s–1。

如图7 所示，为方便查找相同的AIS 报文，将2020 年8 月23 日解析好的AIS 报文导入数据库，选取05:00 UT 时相同的两条AIS 报文，分别为表中第2 行和第4 行。船舶MMIS 码为477108800，经度为128.3945°W，纬度为16.6300°S。

图7 数据库中2020 年8 月23 日天拓五号卫星两个接收机解码的AIS 报文信息Fig.7 AIS message information decoded by the two receivers of Tiantuo V satellite on 23 August 2020 in the database

根据式（8），平面1 与平面3 的夹角为91.0286°，近似垂直。根据式（9），AIS 信号传播的方向向量与平面3 的法向量夹角为28.9486°，近似平行。如图8和图9 所示，为获取功率信息，将解析好的AIS 报文信息导入Excel 文件中显示。如图8 和图9 所示，选取的A 机报文序号为59130，接收功率的16 进制码为02 F3，对应功率值为–111.0814 dBm。选取的B 机报文序号为62641，接收功率的16 进制码为09 E3，对应功率值为–105.8505 dBm。

图8 2020 年8 月23 日天拓五号A 机解码的AIS 报文信息Fig.8 AIS message information decoded by A machine of Tiantuo V on 23 August 2020

图9 2020 年8 月23 日天拓五号B 机解码的AIS 报文信息Fig.9 AIS message information decoded by B machine of Tiantuo V on 23 August 2020

天拓五号卫星上的AIS 天线为±45°交叉极化天线。两天线的全方向图如图10 所示，可以看出两天线方向图对称，且两天线的安装方向与主瓣偏向相反。图10 中坐标轴与图2 中坐标系1 的坐标轴相对应。

根据式（11）可求出1.2节规定的仰角φ为55.0646°。根据式（12）可求出1.2 节规定的方位角θ为90.8541°。对照图10 的天线方向图可以得知两天线的增益近似相等，不用考虑到式（10）的计算中。

图10 双AIS 天线的全方向图（左边为AIS 天线1 方向图，右边为AIS 天线2 方向图）Fig.10 Omni-directional pattern of dual AIS antennas (The left side of the figure is the pattern of AIS Antenna 1,and the right side is the pattern of AIS Antenna 2)

根据式（10）可求出θ1为45.0512°。根据式（14）可求出θ2为28.7046°。根据式（15）可求出法拉第旋转角Ω为73.7559°。根据式（19）可得IPP 坐标为(15.7276°S,123.0625°W)。电离层薄层的高度选取为350 km，通过查询IGRF 可得IPP 处的磁场强度为1.5535×10–5T 根据式（23）可以求出TEC 的值为9.2052 TECU。本次实验使用欧洲定轨中心（CODE）的 TEC 数据与本方法测量的TEC 数据进行对比，CODE 采用全球近200 个GPS 接收站的数据，反演垂直于地表方向的电离层总电子量[18–20]。图11 显示的是2020 年8 月23 日05:00 UT CODE的TEC 全球分布，其中在(15°S,120°W)、(17.5°S,120°W)、(15°S,125°W)、(17.5°S,125°W)位置处的TEC值分别为8.8 TECU、8 TECU、9.7 TECU、8.6 TECU。为使CODE 与本方法的TEC 数据在时间和位置上精确对应，使用双线性插值法对CODE 数据进行时空插值，可以得到05:00:00 UT 时刻CODE 在(15.85°S,123.74°W)处的TEC值为9.1083TECU，对比05:00UT时刻本方法在(15.85°S,123.74°W)处探测到的TEC 值，二者相差0.0969 TECU。

图11 2020 年8 月23 日05:00 UT 的CODE 的TEC 分布Fig.11 TEC distribution of CODE at 05:00 UT on 23 August 2020

CODE 主要是通过陆基观测获取TEC 数据，为确保本方法的TEC 观测环境与CODE 尽可能一致，选取船舶位置位于陆地边缘的AIS 数据进行实验。由于目前该方法还没有实现自动化，找到满足1.1 节条件的AIS 数据耗时巨大，因此本文只找到20 组满足1.1 节条件的AIS 数据进行TEC 测量实验。表1给出了本方法测量的TEC 值与CODE 提供的TEC值的对比。