在轨卫星的工作状态主要通过遥测数据反映, 受工作环境和自身器件退化等因素影响, 地面测控站直接获取的遥测数据中通常含有大量噪声和野值, 因此在分析遥测数据前, 应对其进行滤波处理。目前, 遥测数据滤波方法主要分为基于模型和基于数据驱动两类^[1-2]。当卫星模型已知时, 采用基于模型的滤波方法可以取得很好的滤波效果^[3-5]。但实际卫星系统构成通常较为复杂, 导致精确模型建立困难, 因此, 在工程中多采用基于数据驱动的滤波方法。由于卫星遥测数据包含的噪声和野值在幅值与频率特性上存在差异, 因此在被测数据特性不明确时, 难以使用典型的数据驱动滤波方法同时消除2种干扰。王广成等^[6]和胡婧等^[7]将遥测数据进行窗口化处理，分别采用中值滤波和形态滤波方法，对各窗口内遥测数据进行滤波处理，得到较好的滤波结果，但是结果丢失了数据中的突变特征。周增光等^[8]采用Savitzky-Golay滤波方法对数据中高值噪声和突变干扰具有很好的滤波效果，但计算量较大, 对非平稳数据的滤波能力有限。胡绍林等^[9]提出一种动态测量数据的高保真容错Q算法, 有效消除了遥测数据中的噪声。但方法对野值处理能力有限。可以看出, 现有卫星遥测数据滤波方法对数据应用背景要求较为严格, 在处理噪声和野值问题时, 通常需要分步解决, 增加了计算的复杂度。因此卫星遥测数据的理想滤波是在保留数据中有效数据的同时, 对其中幅值和频率特征不同的噪声与野值进行快速有效滤除。

最速跟踪微分器^[10-11]是由韩京清提出的一种二阶离散形式跟踪微分器, 具有良好的数据跟踪能力, 通过选用合理的参数, 最速跟踪微分器能够实现对含噪数据的理想滤波和去野。但目前最速跟踪微分器的参数选择仍以经验法为主, 并且在滤波时采用全局固定参数, 这都直接影响最速跟踪微分器的滤波效果, 使最速跟踪微分器的应用范围受到局限, 本文提出一种基于滑动最速跟踪微分器的遥测数据滤波方法, 通过滑动窗口对数据进行划分, 根据各窗口内数据的稳定性设计局域最速跟踪微分器的参数, 实现对卫星遥测数据的理想滤波。

1 滑动最速跟踪微分器滤波算法 1.1 滑动最速跟踪微分器参数设计

最速跟踪微分器能够有效抑制输入数据中的高频噪声和野值, 实现对输入数据快速无振荡的跟踪, 因此可以将其应用于遥测数据的滤波处理。本文采用二阶最速跟踪微分器, 其离散形式为

(1)

式中：v(t)为输入数据；x₁(t)为v(t)的跟踪数据和x₂(t)为x₁(t)的导数；fhan(·)为最速控制综合函数；sign(·)为符号函数；r为跟踪因子；h为积分步长；h₀为滤波因子；n为滤波系数, n≥1。

由(1)式可知, 最速跟踪微分器的主要调节参数为跟踪因子和滤波因子。增大滤波因子能够提高滤波输出数据的光滑性, 但也会增大滤波输出数据的相位延迟。文献[12]分析了最速跟踪微分器滤波输出相位延迟与输入数据频率及积分步长间的关系, 证明在输入数据频率确定情况下, 积分步长越小, 则滤波输出数据的相位延迟越小。由于本文滤波对象为离线数据, 因此在实际应用过程中, 通过直接调整积分步长对滤波输出数据的相位延迟, 所以设定滤波系数n=1。为了克服相位延迟, 采用“先预测, 后微分”方法^[11], 对(1)式进行改进

(2)

式中，h₁为预测时长, 其值通常为积分步长h的1~1.5倍。

跟踪因子决定输出数据对输入数据的跟踪程度, 因此通过调节跟踪因子可以直接影响数据滤波效果。如何设定跟踪因子是最速跟踪微分器进行数据滤波关键问题, 目前跟踪因子设定缺乏相关理论, 主要根据经验设定。本文以不同频率和幅值的正弦型数据为输入, 采用扫频法分析最跟踪因子与输入数据之间的幅频关系。

设定最速跟踪微分器滤波因子n=1, 积分步长h=0.01, 输入正弦数据为Asin(ωt), 数据幅值A=[5 10 50], 角频率ω=[5 12 20], 通过组合生成9组不同参数的输入数据。采用扫频法分析跟踪因子与输入数据参数之间关系, 令跟踪因子r初值为1, 步长为1, 分别对各组输入数据进行滤波, 得到输出数据幅值、角频率与跟踪因子之间的关系如图 1所示。

由图 1可得:

1) 随着跟踪因子值增大, 各输出数据幅值均存在转折点, 当超过转折点后, 继续增大跟踪因子值, 输出数据幅值基本保持稳定;

2) 相同角频率下, 不同幅值输入数据对应等幅输出跟踪因子r具有相同的斜率k。

(3)

3) 随着输入数据幅值和角频率参数增大, 对应输出数据幅值转折值逐渐增大, 具体值如表 1所示。

由表 1可知, 相同角频率输出数据的斜率值k分别为:当角频率为5时, 斜率值k_ω5≈0.049;当角频率为12时, 斜率值k_ω12≈0.008;当角频率为20时, 斜率值k_ω20≈0.003。进一步计算各组输入数据角频率与斜率值之间关系可以发现:角频率ω₅=5时, =0.04≈k_ω5; 角频率ω₁₂=12时, =0.007≈k_ω12; 角频率ω₂₀=20时, =0.002 5≈k_ω20, 因此得到斜率k与角频率ω之间存在如下关系

(4)

式中, β为修正系数。由表 1可得到各组数据对应的修正系数分别为:β_ω5=1.225, β_ω12=1.143, β_ω20=1.200, 可以发现β值在一定数值范围内变化, 在实际使用(4)式时, 可设定β=1.200。

综合表 1和图 1结果可知, 对于正弦型数据, 最速跟踪微分器跟踪因子值与输入数据幅值和角频率存在正比关系, 综合(3)和(4)式得到跟踪因子r经验公式为

(5)

目前使用最速跟踪微分器时, 跟踪因子通常设定为全局统一值。当输入数据为非平稳性较强的数据时, 难以在全局内取得较好的滤波结果。为获取理想的滤波结果, 应根据数据的非平稳性征, 对数跟踪因子进行动态更新。

窗口化是处理非平稳数据的常用方法, 通过选择合适的窗口使各窗口内数据的稳定性相对确定。参考(1)式, 为确保最速跟踪微分器相邻数据点滤波结果的连续性, 选择滑动窗方法对数据进行窗口化处理。在选择滑动窗口宽度时, 应使窗口化处理后各窗口内数据呈现相对一致的稳定性。

(5) 式虽然为正弦型数据获得的经验公式, 但其指出跟踪因子值与输入数据幅值和频率之间存在正比关系, 因此参考(5)式, 进一步研究最速跟踪微分器滤波因子与窗口内数据幅值和频率之间关系, 提出一种基于滑动窗最速跟踪微分器的跟踪因子设计方法。

窗口内数据的幅值和频率能够反映数据的稳定性。窗口内数据的幅值稳定性可以利用数据离散程度反映, 而数据离散程度可以通过标准差进行定量表示。窗口内数据的频率稳定程度主要受数据中偏离均值较大的离群点数据影响, 参考广义局域频率定义^[13], 本文提出一种局域最大频率法, 利用各窗口内离群点数据的统计值描述窗口内数据的频率稳定性。方法具体内容为:

首先将原数据划分为M个窗口, 计算各窗口离群点提取阈值

(6)

式中：σ_j(j∈[1, M])为各窗口内数据的标准差；α为阈值系数, 通常取值大于2。

然后统计各窗口内幅值大于离群点提取阈值的数据点数量O_j, 将各离群点数据视为等值的重复数据。最后假设各离群点按照最大频率分布, 可以得到各窗口局域最大频率值为

(7)

式中：N_w为窗口宽度；O_j为窗口内离群点数量。局域最大频率法原理简单, 能够快速计算出各窗口内数据的最大频率。

利用各窗口内数据幅值和频率的定量表征值, 参考(5)式给出各窗口内跟踪因子r_j计算公式

(8)

由于标准差σ_j易受数据幅值影响, 因此采用局域归一化标准差σ_mj对进行修正

(9)

式中, m_j为窗口内的数据均值。将其带入(8)式, 得到各窗口内跟踪因子计算公式

(10)

1.2 滑动窗最速跟踪微分器滤波算法流程

同一卫星系统的遥测数据在不同时段的噪声水平可能存在较大变化。为取得更理想的滤波结果, 在使用最速跟踪微分器进行滤波时, 应对其跟踪因子进行动态修正, 因此本文提出一种基于滑动窗的最速跟踪微分器滤波算法:

首先对输入数据进行滑动窗划分, 滑动窗口宽度设置通常参考输入数据的固有周期, 当数据周期不明确时, 以数据形态特征点, 如拐点划分窗口^[14]。

然后计算各窗口内数据标准差和局域归一化标准差, 设定离群点提取阈值, 统计各窗口离群点数量, 利用(10)式计算各滑动窗内数据对应的跟踪因子, 并使用最速跟踪微分器进行滤波, 得到各滑动窗数据的滤波结果, 以各滑动窗相同数据位滤波结果的均值作为该数据位滤波结果。

由于数据中的突变点可能引起滤波结果产生较大延迟, 因此应对包含突变点的窗口重新划分。数据突变点主要出现在导数值变化较大的位置, 通过设定检测阈值, 提取突变点精确位置, 以此为起始点划分新的窗口, 对窗口内数据的跟踪因子进行更新并重新计算滤波结果, 得到最终滤波结果。方法流程如图 2所示:

2 实验验证及结果分析

为验证滑动最速跟踪微分器的滤波效果, 分别以仿真数据和实测数据进行测试, 并选择滑动中值滤波、滑动多项式滤波和容错Q-滤波进行对比分析。

2.1 滤波效果评价

为了对各种方法滤波效果进行客观评价, 选择信噪比、误差标准差和光滑度作为滤波结果的全局评价参数, 选择突变点平均绝对误差(IMAD)对数据中突变部分的滤波效果进行评价。

信噪比反映滤波器对干扰数据的滤波能力, 其值与滤波器性能成正比; 误差标准差反映滤波结果与输入数据误差的标准差, 其值与滤波器性能成反比; 光滑度反映滤波结果的平滑性, 其值与滤波器性能成反比; 光滑度分为全局光滑度(S)和相对光滑度(S_r), 输入数据{y_i, i∈[1, l]}的全局光滑度计算公式为

(11)

式中, std为标准差计算, Δ_t采样间隔。为滤波结果。

对于包含阶跃跳变的数据, 如矩形数据, 阶跃点会导致数据全局光滑度下降, 因此全局光滑度值无法准确反映滤波结果的光滑度, 针对该类型数据采用相对光滑度

(12)

2.2 仿真数据实验

典型卫星遥测数据的数学形式可简化表示为

(13)

式中：x(t_i)为有效数据；n(t_i)为干扰数据，n(t_i)=n₁(t_i)+n₂(t_i)+n₃(t_i)；n₁(t_i)为随机噪声数据；n₂(t_i)为野值数据；n₃(t_i)为趋势性异常数据。

利用公式(13)生成有效数据x(t)为正弦和矩形的2组仿真数据, 噪声n₁(t)为服从高斯分布, 幅值在[0, 5]区间变化的噪声数据, n₂(t)为数据中均匀分布的10个时长为15的野值。n₃(t)=0.1t为数据中包含的趋势项。

令仿真数据采样时间为0.01 s, 长度N为5 000。当x(t_i)=20sin(5t_i)时, 得到仿真数据y₁(t_i), 如图 3所示。利用本文方法和对比方法对仿真数据y₁(t_i)进行滤波, 统一设定各方法滑动窗宽度为50, 滑动步长为1。数据滤波结果见图 4至7。

各种滤波结果评价参数计算结果如表 2所示。

综合对比图 4至7各方法滤波结果的时域图以及表 2滤波评价参数值, 可知:

滑动多项式滤波结果的误差标准差较大, 全局光滑度和信噪比较低, 对干扰数据滤波能力有限, 滤波结果中仍然包含较多噪声, 突变点平均绝对误差值较小, 对野值几乎没有滤除效果。滑动中值滤波和容错Q-滤波在滤波原理上相似, 2种方法滤波结果的突变点平均绝对误差值较大, 对野值有很好的抑制效果, 滤波结果的全局光滑度接近, 滑动中值滤波具有较低的误差标准差和较高信噪比, 但受到方法原理限制, 滤波结果存在相位延迟和端点效应。与以上各方法对比, 滑动最速跟踪微分器滤波结果的各项评价参数均为最优, 对噪声具有较好的抑制能力。与文献[10]中仿真数据对比, 本文仿真数据包含的野值时长宽度为其3倍, 由滤波结果可知本文方法对宽时长野值仍然具有较好滤除效果。

当时, 令周期T=1 250, 周期数量K∈[0,1,2,3,4], 得到矩形仿真数据y₂(t_i), 如图 8所示。对仿真数据y₂(t_i)进行滤波时, 同样设定各方法的滑动窗宽度为50, 滑动步长为1, 得到滤波结果如图 9至12所示。

各滤波结果评价参数计算结果如表 3所示。

由于仿真数据y₂(t_i)包含周期性阶跃变化, 因此选择相对光滑度作为评价参数, 综合对比图 9至12中各方法滤波结果的时域图以及表 3给出各滤波结果评价参数, 可知:

滑动多项式滤波结果的信噪比和突变点平均绝对误差值较低, 对全局噪声和野值的抑制能力较差。阶跃跳变会引起容错Q-滤波方法中分位数的跳变, 导致对应滤波结果的信噪比和突变点平均绝对误差值下降, 使全局误差标准差增大。滑动中值滤波结果具有较高的信噪比和相对光滑度, 但在阶跃点会产生相位延迟导致滤波结果全局误差标准差增大。阶跃点会引起滑动最速跟踪微分器的局部跟踪参数抖动, 导致相对光滑度下降, 但其他评价参数均为最优。

2.3 实测数据实验

图 13所示为某卫星太阳帆板轴温实测数据, 数据受干扰较为剧烈, 在478点到530点之间出现突变异常。设定各滤波方法滑动窗宽度为30, 滑动步长为1, 各滤波方法的滤波结果如图 14至17所示。

各滤波结果评价参数计算结果如表 4所示。

综合分析滤波结果时域图与相关评价参数, 滑动多项式滤波结果较好地提取了输入数据的趋势, 但对细节保留过多, 导致滤波结果信噪比低、误差标准差较大以及全局光滑度低, 对突变异常数据的滤波能力较差。滑动中值滤波和容错Q-滤波的滤波结果具有较高信噪比, 容错Q-滤波比滑动中值滤波具有更好的光滑度, 但是对突变异常数据的滤波效果有限。与以上方法相比较, 滑动跟踪微分器的滤波结果在全局内误差标准差最小, 具有最高的信噪比和全局光滑度, 能够较好滤除噪声, 提取数据变化的趋势。

3 结论

受工作环境影响, 地面测控站直接获取的在轨卫星遥测数据通常为含有噪声和野值干扰的非平稳数据, 本文提出了一种基于滑动最速跟踪微分器的滤波方法, 利用窗口化处理扩展了最速跟踪微分器在数据滤波领域的工程实用性, 将最速跟踪微分器跟踪因子与数据局部平稳性进行关联, 给出了局部跟踪因子计算方法, 克服最速跟踪微分器使用全局固定跟踪因子的滤波局限性, 能够处理野值和阶跃型的异变数据。采用仿真数据和实测数据进行了实验验证, 结果证明本文方法对不同类型的遥测数据滤波的有效性。