时莎莎，韩昕怡，涂 阔，李一凡，胡劲文

（1.中国空空导弹研究院，洛阳 471009；2.西北工业大学自动化学院，西安 710129）

1 引 言

无人机是发展最成熟、装备量最大、参战最频繁、应用最广泛的无人系统[1]，相对于有人机具有人员伤亡率低、使用条件限制少、隐蔽性好等特点[2]，在民用与军事领域有广泛的用途[3]。尤其是在军事上，无人机与巡航导弹等智能作战设备的结合，无人机的集群作战等是发展的重点方向[4]。但同时，由于无人机上没有飞行员，其起飞和着陆的难度比有人机大得多[5]，因此返航着陆技术成为航空飞行器的关键技术之一，而这一阶段也是在飞机飞行过程中，对导航精度要求最高的阶段。据统计, 世界上飞机进场着陆时的事故发生率约占整个飞行事故的85%，这是因为飞机在进场着陆阶段的飞行比在航线上飞行要困难得多，这一阶段飞机的飞行状态，包括速度、高度、航向等均会发生很大的变化，其中任何一个环节处理不当都有可能引发事故[6-7]。因此，在航空技术的发展过程中，必须对自动着陆系统给予极大的重视。

20 世纪80年代以后，随着技术的发展， 以精确数据为引导信息的“精密进近着陆”引导成为着陆引导系统的发展趋势。仪表着陆引导精度较低, 微波着陆引导系统造价较高，是比较早的着陆系统。同时，地面导航设备也存在着距离导航台越远，误差越大的特点[8]。因此目前发展较为先进的可用于无人机自主着陆的引导技术有雷达着陆引导系统、卫星着陆引导系统、光电着陆引导系统、视觉着陆引导系统等。而这几种自主着陆引导技术也有各自的优缺点。

雷达着陆引导技术具有全天候引导能力，对机载设备要求不高，机动性好，可用于能见度差、云底较低的复杂气象条件下引导无人机着陆。比如美国Sierra Nevada 公司就基于雷达引导技术开发了无人机通用自动回收系统[9]，但是相对比其他系统，雷达引导比较被动，完全依靠地面指挥，获得数据的速率低，引导容量有限[10-13]。光电着陆引导系统使用高分辨率电视、红外热像仪和激光跟踪等复合手段，完成无人机相对于着陆点的偏差数据的测量。具有测量精度高、抗电子干扰能力强，具有高分辨率成像的能力，而且适装性好，很适合与其他引导手段互相匹配，互补共存[14-15]。GPS 系统由于具有传输距离远、覆盖范围广、不受地理条件限制、通信频带宽、能够适应全天候复杂地形等特点，成为目前的主流方式之一。但是卫星通信受自身特点的限制和环境的影响，不可避免地存在各种干扰，特别是其开放式的系统，更容易受到一些不可预见的干扰与专业的软杀伤，比如卫星信号欺骗/干扰等[16]。视觉引导技术是将计算机视觉应用于无人机着陆方式的一种新型引导技术，在卫星通信受到干扰时，可以起到辅助导航作用，已经成为自主导航领域研究的热点[17-18]。视觉导航系统的工作波段远离当前电磁对抗的频率范围，且具有成本低、自主性强、信息量大、无源性好和信息丰富等优点。但是在恶劣天气情况下很难获得比较清晰的图像，从而影响无人机着陆的安全性和精确性[19-22]。

而数据的可靠性不仅仅由导航系统决定，在民航飞机导航系统领域，多传感器耦合导航已成为主流导航方法[23]，其数据融合方式也非常重要。比如早期，Bierman，Carlson 和Schmidt 等于1979年开始先后提出了平方根算法和UD 分解滤波[24-26]，丰富了滤波理论；Seyer 于1979年首次提出了分散滤波的思想[27]；Carlson 于1988年提出了联邦滤波理论，旨在为容错导航提供设计理论[28-29]。其后，各种各样的融合算法不断发展。

由于飞机着陆过程中本身姿态的复杂性，算法的自适应性也非常重要。2008年，Zhao 为了克服噪声统计特性未知时联邦滤波不稳定甚至发散的缺点，提出了一种自适应联邦Kalman 滤波[30]；2009年，Cai 等在解决基于快速块的运动估计中的阈值选取问题时，提出了一种根据目标误差概率的自适应选取阈值的方法[31]；Zhou 等在 SINS/GPS/DVL 组合背景下，采用模糊无迹卡尔曼滤波的联邦滤波算法，进一步提高了联邦滤波的自适应能力和容错能力[32]。随着科技的发展，与信息融合有关的自适应算法不断发展。

基于以上条件，本文针对无人机着陆场景，考虑到数据精度要求高，可用引导源多，不同传感器存在时间延迟等约束，提出了一种无人机着陆过程中的多源引导信息融合算法。该算法重点基于自适应窗长的动态方差估计融合算法，根据导航传感器的量测噪声变化将量测误差方差估计分为平稳段与突变段，通过采取不同窗长的时间窗提高对传感器量测误差方差估计的准确性，并按照误差方差估计值的大小分配融合权重，实现了融合算法的自适应调整，从而提升信息融合精度。

2 数据预处理技术

2.1 时间配准算法

在多传感器数据融合系统的实际工作中，造成数据时间不匹配的原因有以下三点：

（1）多传感器的时间基准一致性问题，即时间同步问题。

（2）各个传感器的采样频率和起始采样时刻不一致的问题。

（3）网络传输延迟造成的时间不匹配问题。在研究中，有许多不同的时间配准方法。比如，Blair等在1991年提出了最小二乘时间配准算法[33]。而本文则选择使用基于内插外推法的时间配准，主要针对量测频率不一致的两组传感器量测信号，以低频信号的量测时间点为参考时序，在同一时间片上对各传感器采集的目标观测数据进行内插、外推处理，将高精度观测时间的数据推算到低精度时间点上，以实现各传感器时间上的匹配[34-37]。

以两组传感器量测数据为例，设传感器A 和传感器B 对同一目标进行观测，在传感器B 的量测时间点Tb0和Tb1间隔内，系统接收到了传感器A 的量测a1X，假设传感器A 在Tai时刻观测到的测量数据为(Xai,Yaj,Zai)，传感器B 在Tbj时刻观测到的测量数据为(Xbj,Ybj,Zbj)，由传感器A 到传感器B 的采样时刻进行时间配准，配准后的数据表示为(Xaibj,Yaibj,Zaibj)

2.2 信号平滑算法研究

信号平滑的主要目的是根据导航传感器误差限，对一定更新范围内统计特性超出误差限的信号进行平滑，防止野点对融合精度产生影响。本节采用基于Savitzky-Golay 滤波的信号平滑方法对超出误差限的各传感器量测数据进行平滑，并作仿真验证。SG 滤波法（Savitzky Golay Filter）的核心思想是对窗口内的数据进行加权滤波，但是其加权权重是对给定的高阶多项式进行最小二乘拟合得到的。它的优点在于在滤波平滑的同时，能够有效地保留信号的变化信息[38-40]。

假设窗长为n=2m+1，同样对当前时刻窗内的2m+1 个观测值进行滤波，各测量点的时间为t=（−m,−m+1,…,−1,0,1,…,m−1,m）用k阶多项式对其进行拟合。拟合公式如下：

在实际验证中取定m=2，k=3，采用五点三次平滑法对超出误差限的量测信号进行平滑。取平滑窗长为5，用于存储量测信息X，最后推算可得

对于初始的前五拍数据，直接进行存储，然后随着时间更新，时间窗也向前移动，每更新一个时间点，平滑窗向前移动一个数据，滑窗内部的最末点为当前时刻的量测信息。

3 多源引导信息自适应融合技术

3.1 技术应用背景

本技术应用背景为无人机着陆时的多传感器融合，因此对无人机着陆情况进行分析。在实际着陆过程中，无人机传感器引导可分为以下四个阶段，如图1 所示。

图1 典型着陆引导系统工作范围示意图Fig.1 Schematic diagram of the working range of a typical landing guidance system

根据着陆过程的不同阶段以及着陆引导传感器的具体工作情况，实际数据融合系统应该能够自动调整融合权重分配，使得融合算法成为一种自适应的融合机制。因此，本文提出了一种基于自适应窗长的动态方差估计融合算法，通过设计自适应时间窗来实时估计传感器的量测误差方差，并根据误差方差的估计结果动态调整融合权重，实现自适应融合。

3.2 基于移动时间窗的残差计算方法

由于实际过程中无法获得无人机的真实位置信息，因此无法获得传感器的量测误差，故本文定义了残差对传感器的量测误差进行估计，并利用实时残差方差对传感器的量测误差方差进行动态估计。在后续过程中，将各传感器残差在特定时间窗内的统计特性关系作为权重信息分配的判断准则，进行数据融合。

3.2.1 残差定义及计算

残差在数理统计中是指实际观察值与估计值(拟合值)之间的差。在本文中，量测值为各传感器实际测得的垂向位置坐标和侧向位置坐标信息与理想下滑道垂向、侧向位置坐标信息之间的差值。由于无法获得真实的传感器量测误差，所以定义残差作为对量测误差的估计，将残差定义为传感器量测值与目标预测值之间的偏差。

令,i kX表示第i种传感器在第k周期的量测值，为第i种传感器在第k周期的纵向偏差量测值,为第i种传感器在第k周期的横向偏差量测值，为第k周期无人机的位置偏差预测值1周期位置偏差的融合估计结果,X5,k与Lk分别是第k周期惯导所测得的位置偏差和理想下滑道的位置，具体计算公式如下：

整理上式后可得

3.2.2 基于拟合度检验的时间窗口选取

由上文可知，残差本身一定符合高斯特性，同时残差数值的大小与时间窗的长度息息相关。时间窗的长度过大，无法保证融合权重分配具有实时性，不一定能够跟上传感器量测噪声的变化。选取的数据量太少，无法体现出所选取时间窗口内的残差数据的统计特性，容易出现过调。如果某几组相邻数据发生跳变干扰，窗长过小会导致融合结果受这几组跳变数据影响增大。

因此，本文选择Jarque-Bera 检验对窗口长度进行判断。具体来说，Jarque-Bera 检验通过计算数据样本的偏度和峰度，评价给定数据服从未知均值和方差正态分布的假设是否成立。

Jarque-Bera 的统计量定义为

其中，n是观测数（或自由度）；S是样本偏度，F是样本峰度。

S，F的计算公式如下：

根据这一特性，我们可以做假设检验。

H0:接受原假设,待检验数据分布服从正态分布；

H1:拒绝原假设,待检验数据分布不服从正态分布。

利用Jarque-Bera 假设检验方法进行多次蒙特卡洛仿真，采用伪随机数的方式进行拟合优度测试，测试不同窗长时间窗内的伪随机数是否符合正态分布，即该组数据是否具有高斯特性从而确定窗口的长度。

3.3 基于自适应时间窗的动态方差估计融合算法

由上文可知，降落过程中飞行状态与传感器精度会产生不断的变化，因此，本文认为融合权重的设置应该能跟随传感器精度的变化而及时做出调整，为了达到这种目的，提出了一种基于自适应窗长的动态方差估计融合算法，具体过程如下：首先，根据传感器量测误差方差的变化幅度将误差方差估计分为平稳段与突变段，从而设置不同的时间窗长，以提高对传感器量测误差方差估计的准确性。接着，利用时间窗内的残差方差作为权重信息分配准则，以实现融合算法的自适应调整。

3.3.1 自适应窗长设计

自适应窗长选取的核心是对传感器量测误差方差变化阶段的判断，即通过检测传感器数据中噪声的变化，将噪声的估计曲线分为平滑段和突变段，对不同的曲线段运用不同长度的窗长对多传感器数据进行融合。

本文采用一维曲线突变段检测方法来检测传感器量测误差方差的变化，由一元连续函数理论可知，信号中的突变段往往对应较大的斜率绝对值[41]。因此可以用导数作为特征,对信号突变段进行检测，如下所示：

总体来讲，自适应窗长设计算法的步骤如下：

（1）采用固定窗n，对各个传感器的量测残差进行方差估计，得到各方差估计点组成的估计曲线。

（2）利用均值滤波对方差信号进行平滑，得到平滑后的方差估计曲线。

（3）利用一维曲线突变段检测算法找出方差估计曲线的平稳段L1和突变段L2。

（4）如果当前判定点x,ik∈L1，把窗长n设置为大窗长；如果判定点属于x,ik∈L2，把窗长n设置为小窗长，窗长由Jarque-Bera 检验确定，对方差进行重新估计，重复步骤(1)、(2)得到新的方差估计曲线。

（5）重复步骤(3)、(4)，根据传感器实时量测误差方差的变化选取不同长度的窗长。

根据以上步骤，可得到不同情况下合适的时间窗窗长。

3.3.2 基于自适应窗长的动态方差估计融合

由前文可得无人机着陆过程可分为四个阶段，针对不同阶段，根据不同的传感器参与情况的量测残差来调整信息融合权重。为了提高精度，将时间窗步长设为s=1，并对算法进行介绍。

着陆第一阶段：无人机进入卫星引导系统的作用范围（d3≤d≤d4）。

由于这一阶段只有卫星引导系统发挥作用，因此经过预处理后的卫星量测即为融合结果：

着陆第二阶段：无人机到达光电、雷达引导系统作用范围（d2≤d≤d3）。

这一阶段，卫星、光电、雷达引导系统均处于有效工作状态，首先根据第k周期（当前时刻t k T= ×Δ）以及第k-1周期的惯导数据X5,k-1,X5,k，结合上一周期的融合结果，按照公式(4)，公式(10)计算第k周期无人机的位置偏差预测值以及三组传感器量测值与预测值的残差,ikr，然后分别计算三组数据在当前时间窗内的残差方差。假设在第1k周期，光电和雷达引导系统开始发挥作用，系统接收到了光电和雷达的量测数据，进入着陆第二阶段，通过以下公式进行融合：

其中，Xk代表第k周期的偏差融合结果，xi,k(i=1,2,3)分别代表预处理后的光电、卫星、雷达传感器的量测值，ri,k代表在第k个周期，第i组传感器窗内残差数据的均值，表达式如下：

Ki,k(i=1,2,3)代表第i组传感器在第k周期时的融合权重，n代表时间窗的长度：

其中，δ为一个远小于1 的小正数，避免方差为0 的情况。代表第k个周期第i组传感器窗内残差数据的方差，表达式如下：

着陆第三阶段：无人机到达视觉引导系统的作用范围（d1≤d≤d2）。

该阶段，四组引导系统都处于有效工作状态，与第二阶段相似，首先按照公式计算与各组残差ri,k(i=1,2,3,4)，然后分别计算4 组数据在当前时间窗内的残差方差，按照以下公式进行融合，r1,k代表在第k周期光电传感器窗内残差数据的均值。

权重表达式如下，Ki,k(i=1,2,3,4)代表第i组传感器在第k周期时的融合权重：

着陆第四阶段：无人机离开光电、雷达引导系统作用范围 （0≤d≤d1）。

此时，无人机到达着陆点附近，离开光电、雷达引导系统的作用范围。这段距离内，只有卫星和机器视觉两种引导方式作用。融合过程表达式如下：

其中，K2,k和K4,k表示经过预处理后的卫星和视觉传感器量测值，K2,k和K4,k代表第i组传感器在第k周期时的融合权重，表达式如下，且满足K2,k+K4,k= 1：

4 仿真验证及结果

4.1 仿真实验场景

为了使实验更加科学准确，根据现实情况，我们设计了如下仿真实验场景。

仿真场景的具体参数设置如下：

（1）无人机匀速向东直线行驶；

（2）无人机东向速度200 km/h、北向速度0、天向速度18 km/h，下滑角4°；

（3）无人机和着陆点的水平距离相距50 km，无人机初始高度4500 m；

（4）量测偏差量（单位：m，最低有效位LSB：0.01 m）；

图2 仿真实验场景Fig.2 Simulation experiment scene

（5）光电、雷达引导手段的有效工作距离为0.1～14.95 km，数据更新频率为20 Hz；

（6）卫星引导手段的有效工作距离为 0～ 50 km，数据更新频率为40 Hz；

（7）视觉引导手段的有效工作距离为 0～ 1.5 km，数据更新频率为40 Hz。

针对上述测试场景，光电和卫星所用的测试数据均来自项目合作单位提供的半物理仿真数据，雷达数据与视觉数据由于传感器厂商无法解算实时数据，所以用于测试的是模拟数据，采用本文所提出的信息融合算法进行仿真测试。

4.2 时间配准算法仿真验证

为了验证时间配准算法，观察配准前后信号的关系，设置两组仿真实验，分别对两组导航传感器采样频率为整数倍以及非整数倍的信号配准进行测试。本文所提到的量测值均为各传感器实际测得的垂向位置坐标和侧向位置坐标与理想下滑道垂向、侧向位置坐标之间的差值。以卫星量测高度信息为例，其更新频率为 40 Hz，即每0.025 s 更新一次，将光电量测的高度信号采用基于内插外推法的配准方法，分别进行整数倍以及非整数倍配准，配准频率分别为 20 Hz 以及30 Hz。

（1）整数倍时间配准

采用内插外推法对卫星高度量测信号进行配准，配准频率定为20 Hz，经过内插外推法时间配准，卫星高度量测的更新频率由原本0.025 s 更新一次变为0.05 s 更新一次。分别计算配准信号与原始信号的误差统计特性得到表1。

表1 时间配准误差统计特性表Table 1 Time registration error statistical characteristic table

（2）非整数倍时间配准

采用内插外推法对卫星高度量测信号进行配准，配准频率定为30 Hz，经过内插外推法时间配准，卫星高度量测的更新频率由原本 0.025 s更新一次变为0.033 s 更新一次。分别计算配准信号与原始信号的误差统计特性，得到表2。

表2 非整数倍时间配准误差统计特性表Table 2 Non-integer time registration error statistical characteristics table

由表1 和表2 可得，经过时间配准后，信号误差均值与误差方差相较于原始信号并无大的波动，配准后的信号能够很好地还原原始信号。相比于整数倍配准，非整数倍时间配准的误差方差的波动幅度结果更大一些，说明非整数配准信号还原原始信号的程度相比于整数倍配准信号的效果略差一些。

4.3 平滑算法模拟仿真

通过仿真验证算法可得，经过数据平滑后，远超出误差限的大幅值干扰信号会被剔除，平滑后的信号噪声明显减小。平滑前后的信号误差统计特性的计算如表3 所示。

表3 光电高度量测信号误差统计特性表Table 3 Statistical characteristic table of photoelectric height measurement signal error

对比平滑前后的信号误差均值与误差方差，可以发现信号误差的统计特性并未发生大幅度变化，能够保留原始信号的噪声特性，具有比较好的平滑效果。

4.4 自适应窗长选取及仿真验证

为验证自适应时间窗设计算法，以光电高度量测信息为例，进行具体讲解。采用一维曲线突变段检测方法来检测传感器量测误差方差的变化，用量测误差方差导数作为特征，对突变段进行检测，根据实际方差估计曲线设置方差导数的阈值检测结果如图3 所示。

图3 光电量测残差方差变化阶段检测图Fig.3 Photoelectricity measurement residual variance change phase detection diagram

图3 中，蓝色曲线为光电高度量测残差方差的导数绝对值，红色水平线为设定的方差导数绝对值的阈值。超出红色水平线的部分判定属于量测噪声突变段，选取较小时间窗长；在阈值水平线之下的采样点判定属于量测噪声变化平稳段，选取较大时间窗长。可以看到光电传感器在初始量测阶段，残差方差变化较大；在采样时间到达875 s 之后，残差方差出现突变，光电传感器在这一阶段的确出现误差变化的突变，说明该检测算法针对这一变化能够进行有效检测。

对于量测噪声变化平稳段和突变段的大小时间窗长的确定，我们利用前文所述的正态分布拟合优度检验方法，对具体的残差时间窗口长度进行选取。

针对残差方差变化的平稳段，结合光电传感器的实际误差统计特性，我们选取方差为0.5 的6 组长度不同的伪随机数序列，分别进行1000 次蒙特卡洛仿真，并对每组的1000 次仿真结果做假设检验，记录其符合正态分布的样本数量。在Jarque-Bera 检验中，通常默认显著性水平为0.05，因此设置显著性水平为0.05。统计结果如表4 所示。

表4 蒙特卡洛仿真结果(平稳段)Table 4 Monte Carlo simulation results (stationary segment)

针对残差方差突变段，结合光电量测实际误差统计特性，我们选取方差为0.25 的10 组长度不同的伪随机数序列，分别进行1000 次蒙特卡洛仿真，统计结果如表5 所示。

表5 蒙特卡洛仿真结果(突变段)Table 5 Monte Carlo simulation results (abrupt section)

据统计结果可以看出，随着时间窗口长度的增加，窗内数据增多，服从正态分布的样本数据也逐渐增多。但根据前文可知，整体上时间窗口的长度不宜设置过大，结合上述统计结果，对于误差方差变化平稳的阶段，我们选择数据长度为200 的时间窗；对于误差方差突变段，本着小窗长的选取原则，选择数据长度100 的时间窗。在实际中，时间窗口的长度为数据长度乘以采样时间。

同样采用上述方法，选定雷达量测误差方差变化平稳段和突变段的时间窗长分别为250 和100。而考虑到卫星传感器鲁棒性较低，与视觉引导系统仅用于近距离着陆，精度要求较高，最终选定卫星量测误差方差变化平稳段和突变段的时间窗长分别为200 和100，视觉时间窗长为100。

4.5 数据融合仿真结果

最后采用基于自适应窗长的动态方差估计融合算法对模拟场景下的飞行航迹进行融合，采用采样频率分别为40 Hz 与20 Hz，由于数据量巨大，为了反映精度，采用分段的方法展示融合结果。按照无人机的飞行高度与不同的着陆阶段进行划分，分别选取纵向高度为1800～2000 m（着陆第一阶段部分航迹）、500～600 m（着陆第二阶段部分航迹）、0～100 m（着陆第三、四阶段部分航迹）的三个高度段，并按照无人机纵向高度对纵向融合误差与侧向融合误差的分布进行绘图，结果如图4 和图5 所示。

由图4 和图5 可得，随着无人机飞行高度的下降，纵向融合误差与侧向融合误差数据明显趋于集中，误差范围也逐渐减小，因此当高度逐渐降低，即无人机逐渐接近着陆点时，会得到更为精确的误差范围，信息融合的精度也会逐渐提高，这符合无人机的着陆需求。

图4 纵向融合误差图Fig.4 Longitudinal fusion error map

图5 侧向融合误差图Fig.5 Lateral fusion error map

5 结 论

本文针对无人机着陆过程中，导航定位精度要求较高，单一导航系统可靠性低、故障率高等问题，提出了基于多传感器的自适应融合导航技术，研究在完成数据初步处理的基础上，根据着陆不同阶段的特点自适应调节算法参数，能够提高融合算法精度，进而提高导航安全性。得到如下结论：

（1）通过内值外插法与基于Savitzky-Golay滤波的信号平滑方法能够有效解决时间配准与野值平滑问题，有效减少信号误差；

（2）通过基于自适应窗长的动态方差估计融合算法，在算法中对各信息源实时残差方差的比值调整融合权重分配，实现融合算法的自适应调整。在无人机接近着陆点过程中，纵向融合误差与侧向融合误差数据明显趋于集中，误差范围也明显减小，有效解决了无人机着陆时可用引导源多、导航精度要求高以及传感器的引导精度时变等问题。