CN109765141B - 一种基于swarm-c卫星提取大气密度的方法 - Google Patents

Abstract

本发明公开了一种基于SWARM‑C卫星提取大气密度的方法，先获取SWARM‑C卫星的加速度计数据和GPS轨道数据，再对其进行重采样和坐标系转换；然后去除加速度计数据中存在的尖峰信号和阶跃信号；数据检查完之后开始对加速度计数据进行校正，用广义最小二乘法估计校正公式中的偏差系数、比例系数和温度系数，获得校正后的加速度计数据；接着对卫星在其运行轨道上受到的光照辐射压进行建模，计算出由于光照辐射压引起的加速度，将校正后的加速度计数据减去光照辐射压引起的加速度，得到大气拖拽加速度；最后计算卫星在其轨道运行的阻尼系数和有效面积的乘积，根据大气拖拽加速度和阻尼系数与有效面积的乘积计算卫星轨道附近的大气密度。

Description

一种基于SWARM-C卫星提取大气密度的方法

技术领域

本发明属于航空航天与空间气候技术领域，更为具体地讲，涉及一种基于SWARM-C卫星提取大气密度的方法。

背景技术

研究大气密度具有重要的意义，它可以应用于不同种类型航天器的精密定轨或宇宙中的其他飞行物体轨迹建模。航天器的精密定轨需要获取其在太空受到的保守力和非保守力，非保守力主要是由大气拖拽力和光照辐射压力组成，而大气拖拽力的大小和大气密度有关。利用大气密度对卫星进行精密定轨及轨道预测可以帮助低轨航天器延长其寿命，也帮助他们重返轨道，分析潜在风险并且在必要时候对卫星发出机动信号。此外，在太空中有大量的在轨飞行物体被归类为空间碎片，因为空间碎片有些具有非常大的面积质量比，若是撞击地球则会造成灾难性的后果，因此需要被仔细跟踪必要时候控制其轨道。譬如在2009年2月两艘完整的航天器发生意外超高速碰撞，该碰撞发生在距地球表面790公里的高度，碰撞后产生的碎片逐渐分离到宇宙中不同的轨道平面，在接下来的几十年中威胁到其他卫星。自从1957年以来，人类已经创造了超过25000个废弃空间碎片，其中许多已经腐蚀，其飞行轨道下降到低层大气中。目前，美国太空监视网络已经追踪了超过20000个大于10厘米的人造物体，那些废弃碎片和卫星碰撞后会使卫星失效，因此对其跟踪和轨迹建模具有非常重要的意义。

大气密度受地磁活动的影响而变化，难以获取其精确值，目前针对热层大气密度的提取方法主要有激光雷达探测法，卫星半长轴衰减法，经典大气模型建模等。大气模型的精度受到地磁活动的影响很大，在地磁活动期间，其精度会下降20％左右。李永平等人在专利《空间用四极滤质器及同时获得大气密度和大气成分的方法》中发明了一种空间用四极滤质器来探测大气密度和大气成分，该方法对仪器制造技术有较高要求，受到硬件技术限制。张鹏宇等人在专利《一种伴飞型大气密度测量装置》中发明了一种包括三向过载传感器、中央处理单元和GPS雷达一体机在内的伴飞行大气密度测量装置。该装置本质还是利用激光雷达，并且其装置的飞行高度有限，无法探测到高层大气特别是航天器所在位置的大气密度。周率等人在《一种热层大气密度修正方法》中通过采用核回归方法计算热层大气密度模型的误差，最后对热层大气密度模型进行修正。然而其修正后的大气模型仍然无法保证地磁活动活跃时期所获得的大气密度结果的精确性。此外还有通过卫星的半长轴衰减法来提取大气密度，该方法通过GPS可以获取到卫星的轨道长半轴数据，然后通过长半轴变化率计算大气密度。然而该方法的动态性不足，其获得的大气密度是卫星某时间段内的平均大气密度。

Swarm系列卫星于2013年11月22日发射，并携带加速度计和GPS接收机作为其科学有效载荷的一部分，通过其携带的加速度计和GPS接收机可以获取大气密度。星载加速度计可以准确而实时地测量卫星所受到的非保守力加速度，卫星受到的非保守力加速度主要由大气拖拽力加速度和光照辐射压加速度组成。因此通过对卫星表面的光照辐射压进行建模获得光照辐射压加速度之后就可以导出大气拖拽力加速度，从而得到大气密度。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于SWARM-C卫星提取大气密度的方法，通过对加速度计数据的准确校正和卫星表面光照辐射压的精确建模，获取卫星轨道周边大气密度的实测数据，提高获得的大气密度结果的实时性和动态性。

为实现上述发明目的，本发明一种基于SWARM-C卫星提取大气密度的方法，其特征在于，包括以下步骤：

(1)、数据重采样

(1.1)、下载SWARM-C卫星的加速度计时间序列和GPS轨道坐标时间序列；

(1.2)、设加速度计时间序列的时间间隔为Δt，GPS轨道坐标时间序列的时间间隔为Δt'；对GPS轨道坐标时间序列进行重采样，使采样后GPS轨道坐标时间序列的时间间隔为Δt，且GPS轨道坐标时间序列和加速度计时间序列长度相同；

(2)、坐标系转换

将地固坐标系下的GPS轨道坐标时间序列转换为惯性坐标系下的GPS轨道坐标时间序列，记为nav时间序列；

(3)、阶跃信号处理

将加速度计时间序列标记为acc时间序列；

将acc时间序列和nav时间序列分别从时刻1开始递增编号，将acc时间序列和nav时间序列在编号为k的位置处的值记为acc[k]和nav[k]；

检测acc时间序列是否存在因为EDAC故障引起的阶跃信号，如果从acc[k]处开始出现阶跃信号，则将acc时间序列中所有t≥k时的值acc[t]减去阶跃的高度h，得到消除阶跃的acc时间序列；

(4)、尖峰信号处理

检测消除阶跃的acc时间序列中是否存在尖峰信号，如果从acc[k]处开始出现了尖峰信号，则剔除acc[k]的值，并在acc[k]处利用插值法重新插值，得到消除尖峰的acc时间序列；

(5)、加速度数据校正

将消除尖峰的acc时间序列记为

在卫星轨道的每个点，定义校准方程为：

其中，

表示nav时间序列数据对时间间隔Δt的二阶微分在xyz三个坐标轴上的投影，B是偏差系数，S是比例系数，Q是温度系数，F是时移常数，G是趋势项系数，t₀是

的起始时间，ε表示噪声相关因子；

通过广义最小二乘法找到校准参数B，S，Q，G的估计值

然后对

进行校正；

其中，

表示对

进行校正后的时间序列；

(6)、计算光照辐射压加速度

使用锥形阴影模型确定卫星在轨道的位置，然后计算出在该位置处卫星所受到的光照辐射压加速度a_solar；

(7)、计算卫星有效面积和阻尼系数的乘积W_k；

设SWARM-C卫星表面共计K块面板，那么，第k个面板的有效面积和阻尼系数的乘积W_k为：

其中，S表示分子热力学最概然分布速度比，γ_k表示第k个面板入射大气分子和卫星面板法向量的夹角的余弦值，A_k表示第k个面板的有效面积，v_out表示大气分子与卫星面板碰撞后的出射速度，v_r为nav时间序列对时间间隔Δt的一阶微分；

(8)、计算大气密度ρ；

其中，m表示卫星的质量。

本发明的发明目的是这样实现的：

本发明一种基于SWARM-C卫星提取大气密度的方法，先获取SWARM-C卫星的加速度计数据和GPS轨道数据，再对其进行重采样和坐标系转换；然后去除加速度计数据中存在的尖峰信号和阶跃信号；数据检查完之后开始对加速度计数据进行校正，用广义最小二乘法估计校正公式中的偏差系数、比例系数和温度系数，获得校正后的加速度计数据；接着对卫星在其运行轨道上受到的光照辐射压进行建模，计算出由于光照辐射压引起的加速度，将校正后的加速度计数据减去光照辐射压引起的加速度，得到大气拖拽加速度；最后计算卫星在其轨道运行的阻尼系数和有效面积的乘积，根据大气拖拽加速度和阻尼系数与有效面积的乘积计算卫星轨道附近的大气密度。

同时，本发明基于SWARM-C卫星提取大气密度的方法还具有以下有益效果：

(1)、本发明针对SWARM-C卫星的硬件特点，增加了尖峰信号处理步骤和阶跃信号处理步骤，剔除了硬件原因引起的加速度计数据误差，使得大气密度提取的结果更加准确。

(2)、本发明针对SWARM-C卫星的加速度数据特点，在传统的加速度计校正公式中添加了温度函数，剔除了温度对加速度计数据产生的影响，使得加速度计的校正结果比传统校正方法准确度提高了15％。

(3)、本发明不受地磁活动情况的影响，即使在磁暴日仍然可以准确地测量热层大气密度，其精度比经典模型提高了20％，与半长轴衰减法相比其动态性更优越。

附图说明

图1是本发明一种基于SWARM-C卫星的大气密度提取方法的流程图。

图2是本发明加速度计阶跃信号处理前后的对比示意图。

图3是本发明加速度计数据校正前和校正后的对比图。

图4是本发明卫星表面光照辐射压建模的时候用到的锥形阴影模型。

图5是本发明卫星各面板阻尼系数和有效面积乘积的时间序列。

图6是本发明得到的大气密度时间分布结果。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

本发明的目的是采用卫星星载加速度计测出卫星在轨道运行所受到的非保守力加速度，再通过非保守力加速度计算大气拖拽加速度，从而获得大气密度。在能够实现的条件下，加速度数据校正尽可能精确，光照辐射压模型建模尽可能完善，这样计算得到的卫星大气密度也会更加精确。为此，本发明提供了一种比较精确的加速度计校正方法，在广义最小二乘法中加入力温度参数的估计，去除了温度引起的加速度计误差，使校正结果精确度得到了提高。

下面我们结合图1，对本发明一种基于SWARM-C卫星的大气密度提取方法进行详细说明，具体包括以下实施步骤：

S1:数据重采样

下载2016年1月日SWARM-C的Level2-daily加速度计时间序列和GPS轨道时间序列数据，加速度计时间序列的时间间隔为30s，对GPS轨道坐标时间序列进行重采样，使得采样后该序列的时间间隔也为30s，重采样后的GPS轨道坐标时间序列和加速度计时间序列长度相同均为540，时间间隔均为30s。此处的加速度计时间序列记为acc。

S2:坐标系转换

将S1中得到的GPS轨道坐标时间序列从地固系转向惯性系坐标。转换公式为：

r_J2000＝Π·Θ·N·P·r_ITRF

其中Π，Θ，N，P分别表示岁差、章动、地球自转和极移对应的旋转矩阵，r_ITRF表示地固系坐标下的轨道坐标时间序列，r_J2000表示惯性系下的轨道坐标时间序列。登录IERS官网查询2016年地球岁差、章动、地球自转和极移对应的旋转矩阵年汇总数据，将S1中的GPS轨道坐标时间序列乘以旋转矩阵得到卫星在惯性系下的轨道坐标时间序列nav。

S3:阶跃信号处理

将S1中的acc时间序列和S2中的nav时间序列分别从1开始递增编号,将acc和nav在编号为k位置的值记为acc[k]和nav[k]。检测acc序列是否存在因为EDAC故障引起的阶跃信号，如果从acc[k]开始出现阶跃信号，那么对acc序列中所有的t≥k时的值acc[t]减去阶跃高度h。阶跃信号校正前如图2(a)所示，阶跃的高度为h，阶跃信号校正后如图2(b)所示，消除阶跃。

S4：尖峰信号处理

检查S3中得到的acc时间序列中是否存在尖峰信号，如果acc[k]出现了尖峰信号，则剔除acc[k]，并且进行插值，使用三次样条插值方法。

S5:加速度数据校正

将S4中得到的未校正的acc时间序列记为

在卫星轨道的每个点，定义校准方程为：

其中，

表示S2中的nav时间序列数据对时间间隔Δt的二阶微分在xyz三个坐标轴上的投影，B是偏差系数，S是比例系数，Q是温度系数，F是时移常数，G是趋势项系数，t₀是加速度序列的起始时间，ε表示噪声相关因子。通过广义最小二乘法来寻找校准参数B，S，Q，G的估计值

接着用下式对曲线计算得到校正后的加速度计数：

其中，

表示校正后acc时间序列。校正后的加速度时间序列如图3。

S6:计算光照辐射压加速度

计算卫星受到的光照辐射压加速度a_solar的时间序列，计算公式如下：

其中，S_⊙表示太阳常数，ξ是比例调整系数，可在卫星定轨时作为参数求解，c表示光速，A_u是一个天文单位长度，r表示卫星和太阳之间的距离，A_k是卫星第k个面板的面积，m表示卫星的质量，n_k表示卫星第k个面板的法向量，s表示卫星指向太阳的单位矢量，ρ_k和δ_k分别表示卫星第k个面板的镜面反射系数及散射系数。f是阴影系数，该系数计算用到了锥形阴影模型，模型中卫星的遮挡物可以使地球或者是月球，该模型如图4所示。如果卫星在本影中f＝0，如果卫星在太阳光照之中f＝1，如果卫星在半影中

其中A表示卫星表面阴影的面积，a表示太阳半径。根据锥形阴影模型计算得到的f的变化范围为0-1，a_solar的量级为10^-9m/s²。

S7:计算卫星有效面积和阻尼系数的乘积

S7.1、计算卫星的单位速度矢量u_D：

其中，v_r表示nav时间序列数据时间间隔Δt的一阶微分，Δt为30s。

S7.2、设SWARM-C卫星表面共计K块面板，计算每个面板入射大气分子流和卫星面板法向量的夹角的余弦值，对于第k个面板余弦值γ_k公式如下：

γ_k＝cos(θ)＝-u_D·n_k

S7.3、计算大气分子热力学最概然分布速度比S：

其中，T代表局地大气温度，m_a大气平均分子质量，这两个参数使用matlab内置标准大气模型计算。R为热力学常数取8314mol^-1K^-1。

S7.4、计算入射大气分子的动力学温度T_in：

S7.5、计算大气分子与卫星面板碰撞后的出射速度v_out：

其中，T_wall表示卫星面板温度，a_E表示入射大气分子与卫星面板的能量交换率。根据ESA公布的SWARM卫星面板材料，此处T_wall和a_E分别取值为270k和0.92。

S7.6、计算卫星的阻尼系数和有效面积的乘积W_k：

其中：

S8:计算大气密度

根据S5中得到的校正后加速度

的时间序列和S6得到的光照辐射压加速度a_solar的时间序列计算大气拖拽力加速度a_atmos的时间序列，公式如下：

然后根据S7中的卫星每个面板的W_k计算大气密度ρ的时间序列，公式如下：

计算得到的卫星阻尼系数与有效面积乘积

如图5，最终得到的大气密度时间分布如图6。

综上所述，本发明所涉及的一种基于SWARM-C卫星的大气密度提取方法，其理论基础为阶跃信号检测与处理，尖峰信号检测与处理，广义最小二乘法参数估计，以及卫星阴影建模。广义最小二乘法参数估计加入了温度参数，去除了温度变化对卫星加速度数据的干扰，比传统的最小二乘法校正结果更加精确，从而提高了大气密度测量的精确度。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于SWARM-C卫星提取大气密度的方法，其特征在于，包括以下步骤：

(1)、数据重采样

(1.1)、下载SWARM-C卫星的加速度计时间序列和GPS轨道坐标时间序列；

(1.2)、设加速度计时间序列的时间间隔为Δt，GPS轨道坐标时间序列的时间间隔为Δt'；对GPS轨道坐标时间序列进行重采样，使采样后GPS轨道坐标时间序列的时间间隔为Δt，且GPS轨道坐标时间序列和加速度计时间序列长度相同；

(2)、坐标系转换

将地固坐标系下的GPS轨道坐标时间序列转换为惯性坐标系下的GPS轨道坐标时间序列，记为nav时间序列；

(3)、阶跃信号处理

将加速度计时间序列标记为acc时间序列；

将acc时间序列和nav时间序列分别从时刻1开始递增编号，将acc时间序列和nav时间序列在编号为k的位置处的值记为acc[k]和nav[k]；

检测acc时间序列是否存在因为EDAC故障引起的阶跃信号，如果从acc[k]处开始出现阶跃信号，则将acc时间序列中所有t≥k时的值acc[t]减去阶跃的高度h，得到消除阶跃的acc时间序列；

(4)、尖峰信号处理

检测消除阶跃的acc时间序列中是否存在尖峰信号，如果从acc[k]处开始出现了尖峰信号，则剔除acc[k]的值，并在acc[k]处利用插值法重新插值，得到消除尖峰的acc时间序列；

(5)、加速度数据校正

将消除尖峰的acc时间序列记为

在卫星轨道的每个点，定义校准方程为：

其中，

表示nav时间序列数据对时间间隔Δt的二阶微分在xyz三个坐标轴上的投影，B是偏差系数，S是比例系数，Q是温度系数，F是时移常数，G是趋势项系数，t₀是

的起始时间，ε表示噪声相关因子，T为acc时间序列的长度；

通过广义最小二乘法找到校准参数B，S，Q，G的估计值

然后对

进行校正；

其中，

表示对

进行校正后的时间序列；

(6)、计算光照辐射压加速度

使用锥形阴影模型确定卫星在轨道的位置，然后计算出在该位置处卫星所受到的光照辐射压加速度a_solar；

(7)、计算卫星有效面积和阻尼系数的乘积W_k；

设SWARM-C卫星表面共计K块面板，那么，第k个面板的有效面积和阻尼系数的乘积W_k为：

其中，S表示分子热力学最概然分布速度比，γ_k表示第k个面板入射大气分子和卫星面板法向量的夹角的余弦值，A_k表示第k个面板的有效面积，v_out表示大气分子与卫星面板碰撞后的出射速度，v_r为nav时间序列对时间间隔Δt的一阶微分；

(8)、计算大气密度ρ；

其中，m表示卫星的质量。

2.根据权利要求1所述的基于SWARM-C卫星提取大气密度的方法，其特征在于，所述的计算光照辐射压加速度的方法为：

其中，S_⊙表示太阳常数，ξ是比例调整系数，c表示光速，A_u是天文单位长度，r表示卫星和太阳之间的距离，A_k是卫星第k个面板的有效面积，m表示卫星的质量，n_k表示卫星第k个面板的法向量，s表示卫星指向太阳的单位矢量，ρ_k和δ_k分别表示卫星第k个面板的镜面反射系数及散射系数，f是锥形阴影模型的阴影系数。

3.根据权利要求2所述的基于SWARM-C卫星提取大气密度的方法，其特征在于，所述的锥形阴影模型的阴影系数的确定方法为：

在锥形阴影模型中选取地球或月球作为卫星的遮挡物，当卫星在本影中时，f＝0；当卫星在太阳光照之中时，f＝1；当卫星在半影中时，

其中A表示卫星表面阴影的面积，a表示太阳半径。

Images