陈唯实 黄毅峰 卢贤锋 张洁 陈小龙
(1. 中国民航科学技术研究院, 北京 100028; 2. 海军航空大学, 烟台 264001)
鸟击航空器是民航的传统安全隐患[1]。 全球每年约发生21 000 起鸟击事件,造成经济损失约12 亿美元[2]。 2009 年1 月15 日发生的全美航空公司1549 号航班哈德逊河紧急迫降事件,是迄今为止最为典型的一起鸟击事故[3]。 随着航班量的持续增长和生态环境的不断改善,中国机场的鸟击防范工作压力持续增大。 2019 年,全国运输机场共发生责任区鸟击事件836 起,同比增加113 起,增幅15.6%,责任区鸟击事件万架次率为0.72,增幅10%。
《运输机场运行安全管理规定》[4]要求机场管理机构应持续开展鸟害防范基础性调研,全面掌握机场内及其附近地区的生态环境、鸟类种群、数量、位置分布及活动规律;绘制鸟类活动平面图;掌握机场内及其附近地区与鸟情动态密切相关的生物类群及影响因素的时间、空间分布情况,分析其中的关系;据此制定和不断完善鸟害防范实施方案,确定各阶段应当重点防范的对象,有针对性地实施鸟害防范措施。
近年来,探鸟雷达技术取得了长足进步,其作为鸟情观测的重要技术手段,能够全天候自动运行并持续积累飞鸟目标雷达数据[5]。 通过对探鸟雷达数据的处理分析,有助于掌握机场周边鸟类活动规律,及时形成鸟情分析报告,指导机场制定科学合理的鸟防措施[6]。 相比之下,传统的机场鸟情生态环境主要依靠人工并借助望远镜和夜视仪完成,调研方法主要包括样带法、样点法和计数鸟群法,且要求调研人员具备鸟类学专业知识[7]。 因此,机场通常会聘请专业的鸟情生态环境调研机构开展相关工作,周期一般为5 年一次。可见,传统的鸟情生态环境调研方法虽然能够识别鸟种,但效率较低、实时性不强,机场管理机构难以及时掌握机场周边最新的鸟情态势,从而制定可实时调整的有针对性的应对策略。 利用先进的探鸟雷达技术,结合人工鸟情调研方法,成为机场鸟情观测的发展趋势。
针对以上需求,本文提出了一种基于探鸟雷达数据的机场周边鸟类目标数量估计方法,采用多目标航迹自动起始跟踪算法对探鸟雷达获取的飞鸟目标信息进行处理分析,实时估计热点区域内的鸟类目标数量,帮助机场管理机构及时掌握机场周边的鸟类数量变化情况及分布态势。
目前,国内外已经研制出了相对成熟的机场探鸟雷达系统。 其中,国外最具代表性的产品包括美国的Merlin 雷达、加拿大的Accipiter 雷达、荷兰的Robin 雷达及英国的Aveillant 雷达,国内的部分高校和科研院所也开展了探鸟雷达技术研究,探鸟效果逐步提升[8]。
图1 为典型的双天线机场探鸟雷达系统(https://detect-inc. com/)。 其通常配有2 部不同波段的雷达,一部为X 波段垂直扫描雷达,另一部为S 波段水平扫描雷达,均采用T 型波导阵列天线,大部分早期成熟的探鸟雷达均采用这种体制,如Merlin、Accipiter 和Robin。 水平扫描雷达能覆盖机场周边的广大地区,并获取水平平面的二维信息,但不能提供高度信息;垂直扫描雷达可以获得垂直平面的高度信息。 由于水平和垂直扫描的2 部雷达覆盖的范围不同,波束重叠区域很小,系统获得的三维信息不是真正的三维。 探鸟雷达系统通常被放置在机场跑道的附近或一端,对包括航班起降通道在内的机场区域进行监控。 垂直扫描雷达与跑道连成一线,使其扫描范围恰好覆盖航班的起降通道,同时提供该区域的鸟类目标高度信息;水平扫描雷达对机场周边区域进行扫描,对进入该区域的鸟类目标进行预警。 对于小型鸟类,该系统的探测范围可以达到垂直高度2 300 m,水平距离2 ~3 n mile(1 n mile =1.852 km);对于中大型鸟类或鸟群,探测范围可以达到垂直高度15 000 ft(1 ft =0.304 8 m),水平距离为4 ~6 n mile。近年来,以Aveillant 雷达系统为代表的相控阵、全息雷达等先进的雷达技术逐步应用于机场探鸟,进一步提高了鸟情数据的可靠性,引领了探鸟雷达技术的发展方向。
图1 典型的双天线机场探鸟雷达系统Fig.1 Typical airport avian radar system with dual antennas
雷达获得的量测数据包括目标信息及部分杂波干扰,跟踪目的在于从杂波中提取出每个目标形成的航迹。 目标跟踪要先解决目标状态的初始化关联,即航迹起始问题[9]。 现有的航迹起始跟踪算法主要包括顺序处理和批处理两大类[10]。顺序处理算法通常基于一定的逻辑规则,主要包括直观法[11]和逻辑法[12],此类算法计算量较小,但在杂波环境中性能下降明显。 批处理算法主要包括Hough 变换法及其衍生的改进算法[13-14],此类算法发源于图像处理,可以起始直线运动的目标且适用于杂波环境,但由于计算量较大导致航迹起始的实时性欠佳。 近年来,支持向量机[15]、随机森林[16]、神经网络[17]等智能算法也被用于航迹起始,但此类算法通常只适用于特定的训练数据,普遍存在可移植性差和时效性不佳的问题。
针对以上问题,本文结合机场鸟类活动规律,基于探鸟雷达获取的鸟情数据,以鸟类目标活动热点为参照点,建立飞鸟目标航迹起始与消亡的概率分布模型,缩短了航迹起始时间,降低了航迹虚警率。 本节首先介绍多目标航迹自动起始跟踪的算法流程,进而讨论多目标数据关联算法、目标生命周期中各类可能事件的关联概率估计及卡尔曼滤波与平滑方法。
多目标航迹自动起始跟踪算法的关键和难点在于目标状态的初始化。 基于当前时刻的状态预估值,由粒子滤波的数据关联算法对当前时刻的量测值与可能事件一一关联,包括目标的新生、延续和消亡,同时排除杂波影响。 在量测值与所有已知目标的预估状态和新生目标的初始状态关联之后,进行卡尔曼状态更新,获得当前存在的所有目标的状态更新值。 对全部卡尔曼滤波结果进行平滑处理,得到每个存在过的目标的平滑航迹。算法流程如图2 所示。
图2 多目标航迹自动起始跟踪算法流程Fig.2 Flowchart of multi-target automatic initiation and tracking algorithm
每个目标都有其生命周期,即从新生到消亡的过程,因此目标存在的数目始终处于变化之中。通过数据关联算法,判断量测属于已存在的目标、新生目标、抑或杂波,从而将复杂的雷达多目标跟踪问题简化为普通的单目标跟踪问题,同时记录每个目标的生命周期,对当前存在的目标总数进行实时统计。
基于粒子滤波算法进行的多目标数据关联,可视为多假设跟踪(Multiple Hypothesis Tracking,MHT)算法的推广。 每个粒子代表不同的数据关联假设,降低了MHT 算法的复杂性。 在每个粒子中,将量测与若干可能事件相关联,计算每个事件的先验概率Pp和事件关联相似度Pl,进而给出每个事件的关联概率P,2.3 节将列出所有可能事件的关联概率模型。
估计有效样本数目ne,且
如果ne过小,则进行重采样。
在没有任何先验知识的情况下,量测的关联结果需要基于对杂波、新生目标、已知目标、目标消亡等所有可能事件出现概率的计算。
设pb为航迹起始即新生目标概率,新目标的出现意味着当前量测作为新目标的起始状态。 本文中,将飞鸟目标活动热点作为目标航迹起始的参考点,通过量测与热点的关联度起始航迹,降低了跟踪算法的复杂性,并提高了其效率与工程适用性。 另外,设pd为目标消亡概率,以gamma 分布建模:
式中:θ为目标航迹与量测关联的中断时间,中断时间越长,目标消亡的概率越大;函数的坡度随α和β的不同取值而变化。
因此,建立所有可能的量测关联事件概率模型,并归纳如下:
1) 杂波,无目标消亡。
5) 新目标航迹起始,无目标消亡。
基于以上6 类可能事件的关联概率估计结果,实现对每个目标从新生、延续直至消亡的全生命周期管理。
卡尔曼滤波的实质在于:根据目标的当前状态估计其未来的状态,并通过量测值对其进行修正。 卡尔曼滤波分为2 步:①预估。 根据前一步的量测值预估出系统下一步的状态。 ②更新。 根据当前时刻的量测值估计出系统的当前状态。 其中,预估部分由方程(13)完成:
式中:mk和Pk分别为k时刻获得量测之后估计的状态均值和方差;v为量测修正;K为滤波增益,其定义了预估值修正的程度。
卡尔曼平滑与卡尔曼滤波的流程类似,区别在于滤波的递推是前向的,而平滑的递推是后向的。 经过卡尔曼平滑的处理,可以获得比仅进行卡尔曼滤波更高的精度。
本节针对不同杂波环境中的目标跟踪仿真数据,通过蒙特卡罗实验,评价本文算法特别是多目标航迹自动起始跟踪算法的有效性,重点在于不同参数设置情况下目标航迹起始的及时性及目标数目的变化情况。
运动目标状态表示为
式中:(xk,yk)为二维直角坐标系中目标在k时刻的位置;(̇xk,̇yk)为速度。
目标的离散动态模型可以表示为线性、时不变构造方程:
其中:时间步长设定为Δt=0.01,过程噪声的功率谱密度设定为q=0.001。
图3 为杂波环境中目标的仿真航迹,仿真步数为50 步。 量测中属于目标的部分附加了一定的高斯噪声,杂波均匀地分布在[ -5,5] ×[ -1,9]的区域内,其在每个扫瞄周期中出现的数量满足均值为λ的泊松分布。 以(0,0)为目标起始的参考点,目标航迹的起始点与参考点的距离为d,目标起始速度的绝对值为1,起始方向在0° ~360°范围内随机分布。 图3 为某次仿真实验的目标航迹示意图,图3(a)为全局示意图,图3(b)为局部放大图,该实验中d=0.2,λ=0.5,且粒子数N=50。
图4 为另一次仿真实验的目标航迹示意图,图4(a)为全局示意图,图4(b)为局部放大图,该实验中d= 3,λ= 1,且粒子数N= 50。 对比图3和图4 的2 次目标跟踪结果可知,图3 中的目标航迹起点距参考点很近(d=0. 2),航迹起始没有延时;图4 中的目标航迹起点距参考点相对较远(d=3),航迹起始出现了1 个周期的延迟。
图3 无延迟起始的目标跟踪起始仿真数据Fig.3 Simulation data of target tracking initiation without initiation delay
图4 延迟起始的目标跟踪起始仿真数据Fig.4 Simulation data of target tracking initiation with initiation delay
通过设置不同的杂波环境(λ值)和目标航迹起始条件(d值),图5 和图6 统计了100 次蒙特卡罗仿真实验的目标跟踪结果的平均值,考察了不同参数设置条件下目标跟踪的稳定性与航迹起始的及时性。
图5 给出了目标航迹起始位置与参考点不同距离情况下(d=0.1, 0.5, 1.5, 2.5, 3)目标航迹起始延迟周期随杂波环境的变化情况。 可见,随着杂波数量的增加(λ=0 ~1.5),目标航迹起始的延迟周期波动性增加,较多的杂波影响了目标关联的及时性,导致了航迹起始的延迟。 延迟周期曲线的波动性可能与蒙特卡罗实验的次数有关,其波动性会随着实验次数的增加而降低。由图5 的实验结果可知,当d<3(d=0.1, 0.5,1.5, 2.5)时,航迹起始延迟周期曲线的取值范围变化不大,总体上小于1,明显优于3 个扫描周期完成航迹起始的逻辑法[9]。 即使当d= 3时,航迹起始延迟周期曲线的取值仍然总体小于3。
图6 给出了不同杂波环境下(λ=0.1, 0.5,1.0, 1.5, 2.0)目标航迹起始的延迟周期随航迹起始点与参考点距离的变化情况。 与图5 所得结论类似,随着相对距离的增加(d=0 ~2.5),每条航迹起始延迟周期曲线具有一定起伏,但总体取值范围的变化不大。 另外,当杂波较少时(λ=0.1),航迹起始延迟周期接近于零,几乎没有延迟;当λ取值逐渐增大时,航迹起始延迟周期随之增大,但总体不大于3。
图5 不同杂波条件下的目标跟踪延迟周期Fig.5 Target tracking delay period under different clutter conditions
图6 不同起始条件下的目标跟踪延迟周期Fig.6 Target tracking delay period under different initiation conditions
表1 给出了λ不同取值情况下航迹起始延迟周期曲线的平均值,全部小于2,低于3 个扫描周期完成航迹起始的逻辑法。 当λ=2.0 时,本文算法所得的航迹起始延迟周期平均值为1.858 3,这是由仿真数据中混入的杂波造成的,本文建立的仿真模型中用一部分杂波数据替换了目标量测数据,导致了目标航迹一定程度的延迟起始。
表1 目标航迹起始延迟周期平均值对比Table 1 Comparison of delay periods of target path initiation
如图7 所示,本实验设计的仿真总步数为500,时间步长为0.01,对依次起始的3 个目标航迹进行了仿真。 目标1 的生命周期为t=0 ~5,其从( -0.25, -1.5)出发,在t=0 ~2 以速度(0,1)匀速运动,在t=2 ~4 完成左转弯,在t=4 ~5以速度( -1,0)匀速运动,直至航迹结束;目标2的生命周期为t=1 ~4.5,其从(0, -1.5)出发,在整个生命周期中以速度(0,1)匀速运动,直至航迹结束;目标3 的生命周期为t=0.5 ~4.5,其从(0.25, -1.5)出发,在t=0.5 ~3 以速度(0,1)匀速运动,在t=3 ~4 完成右转弯,在t=4 ~4.5以速度(1,0)匀速运动,直至航迹结束。 目标的真实运动轨迹如图7 中实线所示,“ 。”代表目标量测,其附加了一定的高斯噪声,“ ×”代表杂波,其参数设置为λ=0.1,均匀地分布在[ -3,2] ×[ -2,3]的区域内。 在数据关联中,粒子数设定为N=50,pb设定为0.01,cp设定为0.02,cd设定为1/16,目标航迹的起始参考点设置为(0,0)。
图7 和图8 对比了在目标消亡模型设置不同参数的情况下3 个目标航迹的跟踪结果和目标数目估计结果,2 组实验都实现了杂波环境中全部目标航迹的及时起始和稳定跟踪。 在航迹终结时间估计方面,当目标消亡模型中的参数设置为α=2 和β=0.5 时,其对目标航迹起始与消亡时间的判断与真实情况基本相符,如图7(a)和图8(a)所示。 当目标消亡模型中的参数设置为α=2 和β=2.5 时,其对目标消亡时间的判断出现延迟,目标2 和目标3 的航迹在t=4.5 均未能及时结束,导致目标数据的估计结果出现偏差,如图7(b)和图8(b)所示。 可见,参数β通过改变目标航迹与量测关联的中断时间控制目标消亡概率,β值越大,目标消亡概率越低。
图7 多目标跟踪仿真实验Fig.7 Multi-target tracking simulation experiment
图8 目标数目估计Fig.8 Target number estimation
本节基于探鸟雷达在国内某机场获取的鸟情数据,采用多目标航迹自动起始跟踪算法,对机场周边的鸟类目标数量进行统计分析。
图9 为中国民航科学技术研究院探鸟雷达实验系统。 该系统对S 波段导航雷达进行了升级改造,采用可靠性更高且发射功率较小的固态发射机替代原有的磁控管,引入脉冲压缩和脉冲多普勒等多项信号处理技术,开发了专门的鸟情数据处理与统计分析软件,改善了其在复杂低空环境中对飞鸟目标的探测能力,该系统放置于某机场灯光站附近,已获取并积累了大量机场周边鸟情信息。
图9 中国民航科学技术研究院探鸟雷达实验系统Fig.9 CAST experimental avian radar system
机场周边的鸟类活动通常遵循一定规律,利用探鸟雷达获取的机场及周边区域鸟情信息,对机场留鸟的活动节律进行统计分析,能够掌握机场周边的鸟类数量与分布情况。 本节基于中国民航科学技术研究院探鸟雷达实验系统获取的国内某机场飞鸟目标数据,给出了机场周边鸟类目标数量进行统计分析的基本步骤及实例。
1) 鸟类活动时段统计
基于探鸟雷达获取的该机场某日的鸟情信息,统计探鸟雷达监视范围内全天中每小时出现的飞鸟目标量测数量,经统计,第8 个时段(7:00—8:00)和第18 个时段(17:00—18:00)出现了2 次鸟类活动高峰。 结合机场周边的鸟情人工调研结果,该机场周边的鸟类通常会在机场周边筑巢,并在特定时间进入飞行区觅食。 以上雷达观测结果与人工调研结果吻合,这2 个时间段可以判定为周边鸟类的离巢和归巢时间。
2) 鸟类活动网格分布统计
在探鸟雷达的监控范围内划分网格,该区域面积为L×L,单位为m2。 每个网格的面积为l×l,单位为m2。 本例中,探鸟雷达的覆盖半径为3 km,L=6 000 m,l=100 m,监控范围内网格总数为60 ×60,如图10 所示,左上角设置为坐标系原点(0,0),X轴水平向右,Y轴垂直向下。 在机场周边的鸟类离巢时间段(7:00—8:00)内,鸟类目标量测超过一定阈值的网格均做颜色标记。 对于标记后的网格分布图进行腐蚀和膨胀等二值图像处理,消除其中的孤立点,形成图10 所示的连通区域即热点,结合机场周边的鸟情调研结果,可判定为鸟类的栖息地。 以该热点区域的中心点作为目标起始的参考点,本例中,该参考点坐标为(4 750,2 650) m。
图10 机场周边鸟类活动热点统计结果Fig.10 Hot spot statistics of bird activities around airport
3) 鸟类数量统计分析
图11 为某日清晨鸟类飞离热点区域内栖息地的飞行轨迹分布情况,包括雷达覆盖范围内的全局图和鸟类飞行轨迹的局部放大图。 该鸟类栖息地位于飞行区外东北方向,图11 的离巢过程发生在7:00—8:00 之间,持续约1 min,种群数量在20 只左右,飞行方向大体为西南,通常会进入飞行区取食。
图11 鸟类飞行轨迹分布Fig.11 Flight path distribution of birds
图12 为不同日期清晨7:00—8:00 对该栖息地内离巢鸟类目标数量进行估计的结果,目标跟踪时间持续约70 s。 在数据关联中,粒子数设定为N=50,pb设定为0.01,cp设定为0.02,cd设定为1/16,目标航迹的起始参考点设置为(4 750,2 650)m。 目标消亡模型中的参数设置为α= 2和β=0.5。 以每次鸟类目标数目估计过程中的最大值作为该种群数量的估计值,因此,第1 次估计结果为23 只,第2 次估计结果为22 只。 表2给出了该栖息地连续10 天的鸟类目标数量统计结果。 可见,该鸟类的种群数量约20 ~30 只,且与人工观测结果基本吻合。
图12 鸟类数量估计Fig.12 Bird number estimation
表2 鸟类目标数量估计结果Table 2 Estimation results of bird target number
本文基于粒子滤波数据关联算法,通过建立目标的新生与消亡模型,较好地实现了杂波环境中雷达多目标的自动起始跟踪,并将其成功应用于机场周边的鸟类目标数量估计问题,得出以下结论:
1) 本文提出的自动起始跟踪算法在数据关联之前设定了一个公共的新生目标起始关联点,通过计算量测值与该起始点的关联概率,判断其是否属于新生目标。 该起始关联点的选择基于对目标运动轨迹起始范围的掌握,在一定程度上依赖先验知识。
2) 通过仿真实验证明,本文算法在目标起始的及时性方面明显优于3 个或4 个扫描周期起始的逻辑法,即使在杂波环境中,其目标起始的延迟周期仍然小于2。
3) 本文算法的局限性在于其适用于目标起始区域范围相对集中的多目标自动起始跟踪问题,对于目标起始点分散的情况,工程应用中仍推荐3 个或4 个扫描周期起始的逻辑法。
4) 将本文算法应用于机场周边的鸟类数量统计分析,能够掌握机场周边鸟类的栖息地分布、种群数量及活动规律,进而指导机场开展有针对性的生态环境治理措施,有效提升机场鸟击防范科学水平。