气象雷达探鸟技术综述

2023-04-19 04:30陈唯实刘佳王青斌卢贤锋张洁陈小龙黄毅峰
航空学报 2023年5期
关键词:鸟群鸟类气象

陈唯实,刘佳,王青斌,卢贤锋,张洁,陈小龙,黄毅峰

1.中国民航科学技术研究院 机场研究所,北京 100028

2.北京航空航天大学 前沿科学技术创新研究院,北京 100191

3.海军航空大学,烟台 264001

随着中国民航运输的快速发展和生态环境的持续改善,鸟类的生存活动范围不断扩大,鸟击作为影响航班运行安全的重大隐患之一,对飞行安全的影响持续增加[1]。据中国民用航空局统计数据,中国民航鸟击不安全事件依然呈高发态势,严峻的鸟击防范形势要求机场结合鸟情发生规律有针对性地开展防治工作。春秋季为中国民航鸟击事件高发期,季节规律明显,候鸟迁徙季节中发生在中高空及夜间的鸟击尤为难防难控,其发生规律与候鸟迁徙节律高度相符[2]。

在全世界范围内,每年有数以亿计的鸟类通过长距离迁徙满足其生存和繁殖要求,包括迁移到季节性最佳的栖息地、觅食、躲避捕食者、寻找资源和配偶等。研究鸟类迁徙极具挑战性,因为迁徙鸟类通常会超出观测者的视距范围,必须依靠技术手段进行跟踪[3]。卫星跟踪是鸟类迁徙研究的重要手段,卫星定位设备的小型化和技术革新使其跟踪能力持续提升,然而这些设备的相对重量仍然是限制该技术路线发展的主要因素[4];因为只有一小部分体型较大的鸟类可以通过安装卫星定位传感器进行跟踪,且通常每次被定位的个体数量极为有限,只占种群总数中很小的一部分。雷达便于量化生物目标在较大空间范围内的分布和运动情况,20 世纪40年代已被应用于鸟类迁徙研究[5]。近年来,随着雷达技术的快速发展,已经出现了专业的探鸟雷达,但此类雷达探测范围一般较小且部署数量有限,难以通过组网实现对鸟类迁徙的观测。

气象雷达特别适于大尺度范围内的迁徙鸟群目标探测,且便于有效组网[6]。因此,通过建立基于气象雷达的全国鸟情预警系统对候鸟迁徙情况进行监视预警,将有助于降低目前居高不下的鸟击事件数量。本文在对现有世界范围内的气象雷达组网鸟情探测预警技术进行介绍的基础上,分析了气象雷达获取的飞鸟目标回波特征和气象雷达鸟情信息提取技术,进而讨论了基于气象雷达的鸟类学研究与应用情况,最后分析了在中国建立基于气象雷达组网的全国鸟情预警系统的可行性,并做出结论。

1 鸟情监视预警系统

雷达数据能够提供飞行中生物目标的空间分布、密度、高度、速度和方向等信息,双极化雷达甚至能提供目标的形状信息,极大提高对鸟类、昆虫、蝙蝠和降水等目标的识别能力[7]。自20 世纪40年代起,基于雷达的鸟情监视预警技术已经从零散的研究发展为有组织的研究。目前,美国的下一代气象雷达(NEXt generation weather RADar, NEXRAD)[8]和欧洲的气象雷达信息交换运行计划(Operational Program on the Exchange of weather RAdar information, OPERA)[9],已经能够通过气象雷达组网的方式获得大陆尺度上的鸟类迁徙数据,以建立起气象学家与生物学家之间的有效合作渠道,为鸟类生态学研究提供了独特而广泛的信息来源,有助于更好地理解鸟类迁徙行为,为相关生态学研究以及鸟击防范工作创造条件。

1.1 早期预警系统

20 世纪中叶开始,雷达就被用于鸟类迁徙行为的记录和监视。雷达的主要优势在于,其能够不受光照条件影响,实现对飞鸟目标的全天候监视,记录高空鸟类的数量,并测量其飞行方向和速度,但缺点是不能识别鸟种[10-11]。最早的鸟情预警系统应用于军机保障,其基于对空监视雷达的显示器图像提取鸟情信息。在荷兰,雷达站点会委托生物学家确定雷达屏幕上出现的鸟群密度[12]。在丹麦、比利时和德国,雷达操作者通常自行分析雷达图像。在以色列,这一任务在鸟类迁徙季节由野生动物管理者和雷达专业技术人员共同完成,一旦发现大量鸟群,会立即通知空管塔台及相关地区航路上的空管单位[13]。

最早的鸟情预警电子系统于1971年在丹麦和1978年在荷兰投入使用[14],这2 个系统都是军用对空监视雷达屏幕上多个窗口嵌入的鸟类迁徙强度电子计数器。该类系统进行技术升级的关键是研发鸟类目标的检测、跟踪与识别算法,进而实现飞鸟目标与飞行器目标的分类,并通过定制的软件展示出飞鸟目标的速度、高度和飞行方向。此类升级后的鸟情预警系统已于20 世纪90年代在德国和荷兰部署。其中,德国的系统已经发展成一个由19 部空管雷达和32 部来自德国及周边国家的对空监视雷达组成的网络[15]。21世纪之后,专业探鸟雷达逐渐出现,荷兰的罗宾探鸟雷达能够基于军用对空监视雷达实现鸟情探测,区分出飞鸟轨迹、地杂波和降水信息[16]。荷兰代尔夫特理工大学和德国航空航天中心的研究人员利用探鸟雷达获取的鸟情数据,提出了鸟击风险评估算法,为机场鸟击防范工作提供了有益参考[17]。表1 梳理了鸟情监视预警系统的早期发展情况。同其他雷达一样,气象雷达也能够实现对飞鸟和昆虫等生物目标的探测,尤其便于组网以实现较大尺度上迁徙行为的监视[18]。目前较为成熟的气象雷达鸟情预警网络包括美国的NEXRAD 系统和欧洲的OPERA 系统。

表1 鸟情监视预警系统的早期发展情况Table 1 Early development of bird monitoring and surveillance system

在迁徙预警模型方面,大量研究表明,鸟类的迁徙行为与气象因素有关[19]。迁徙行为与气象条件的关系促进了鸟类迁徙密度与飞行高度预估模型的建立,根据气象预报数据,可以对鸟类迁徙进行24~72 h 的预测。20 世纪70年代初,丹麦率先开发了飞行安全预测模型,但未投入实际应用[20]。德国开发了首个实际应用的飞行安全模型,该模型使用多个天气参数预测鸟击风险。目前,该预测模型还结合了雷达系统的实际鸟类迁徙信息,对预测结果进行了微调[21]。根据气象数据和部署在荷兰和比利时的2 台军用对空监视雷达对迁移强度的测量数据,荷兰建立了一个回归模型,并于2001年投入运行[22]。随后,开发了整体预报模型并投入使用,每小时预报每个雷达站的迁徙情况,该首次尝试仅基于一年的雷达数据,在一个框架内使用气象雷达来预测多个地点的迁移强度及飞行高度分布,因此不足以满足实际应用需求。

1.2 美国NEXRAD 系统

美国的全国鸟情预警系统(Avian Hazard Advisory System, AHAS)由DeTect 公司为美国空军研发并运营,系统能够自动处理来自美国国家气象雷达网络中的NEXRAD 数据,并为美国(包括美国大陆、阿拉斯加、夏威夷、关岛和韩国)所有的空军训练航线及活动区域提供实时航空器鸟击风险评估信息[23]。NEXRAD 气象雷达网络在美国境内的160 个S 波段WSR-88D 多普勒天气雷达站点位置见图1[8],该气象雷达网采集的信号由DeTect 公司位于佛罗里达州巴拿马城的总部通过卫星接收、并经过汇总分析处理后识别出鸟情态势信息并进行发布。

图1 NEXRAD 气象雷达网点分布及覆盖范围[8]Fig.1 Weather radar distribution and coverage of NEXRAD[8]

自1957年美国气象监视雷达基础设施(WSR-57)初步部署以来,对雷达网络的几次连续升级在增强其气象观测能力的同时也为鸟情预警提供了更好的数据基础[24]。早期的气象雷达仅能提供雷达反射率因子,其在1988年升级到WSR-88D 雷达之后才具备了多普勒探测能力,由此获取的多普勒径向速度和径向速度谱宽2 个参数,通常被用于描述迁徙目标的回波特征、区分鸟类和昆虫目标、以及识别鸟类栖息地的回波[25]。2013年,NEXRAD 气象雷达网完成了全部WSR-88D 雷达的双极化技术升级,增加了总差分相位、差分反射率和共极相关系数3 种常规数据产品,并建立了含有目标极化特征的标准数据格式,相关数据产品可通过亚马逊云计算服务平台免费下载,为迁徙鸟类探测与相关研究创造了优越条件[26]。

1.3 欧洲OPERA 系统

在欧洲,所有运行中的气象雷达持续监测大气状况并在OPERA 系统中进行组网,以便在欧洲大陆范围内进行数据交换和测量,该系统于1999年建立,已集合了30 个成员国的200 部雷达,如图2 所示[9,27]。从2013—2017年,欧洲动物活动雷达监视网(European Network for the Radar surveillance of Animal Movement, ENRAM)在欧洲科技合作计划的资助下建立并在OPERA 系统的组织下运行,其主要目标是融合专业知识,在欧洲建立一个包括生态学、鸟类学、昆虫学、气象学、工程学、数学、视觉分析与计算机科学在内的多学科交叉的研究网络,利用OPERA 天气雷达监测动物在欧洲各地的迁飞情况,并提供相关服务。目前,ENRAM 已经采用了OPERA 系统中覆盖西欧部分地区的70 部气象雷达,且这一数字仍在持续增长[28]。

图2 OPERA 组网雷达的地理分布[9,27]Fig.2 Radar distribution and coverage of OPERA[9,27]

欧洲的雷达网由各国的国家气象部门独立建立,并于20 世纪90年代开始了经验交流、协调和数据交换工作[29]。欧洲各国处于其雷达运行网络发展的不同阶段,在雷达翻新或更新方面有不同的战略,半数以上雷达建于1998年之前,网络中最老和最新雷达之间的时间跨度超过40年。在OPERA 系统目前的200 部雷达中,184 部具有多普勒功能,106 部为双极化,且这一数字仍在稳定增长。另外,这些雷达中的大部分工作在C 波段,只有欧洲南部部署的一部分雷达为S 波段。因此,欧洲的OPERA 雷达网络在安装日期、制造商、扫描策略、信号处理和产品输出等方面呈多样化,与美国的NEXRAD 雷达网截然不同,表2对二者的主要参数做了对比[9]。这对于一个大陆网络的运行来说,无疑增加了困难和挑战,要求国家间开展合作与信息交流,因此,OPERA 系统已于2011年成立了称为“奥德赛”的OPERA 数据中心,并建立了数据信息模型(OPERA Data Information Model, ODIM),但由于各国的雷达网仍然独立运行,这种差异短期内难以完全弥合[30]。

表2 NEXTRAD 与OPERA 网络的主要参数对比[9]Table 2 Comparison of main parameters between NEXTRAD and OPERA network[9]

ENRAM 项目重要的功能之一,是通过飞行安全鸟类规避模型(FlySafe Bird Avoidance Model, FlySafe-BAM)提取鸟情预警信息并发布[31]。FlySafe-BAM 由欧洲的多家科研机构合作研发,实现了对鸟类迁徙的自动检测和量化,能够给出不同高度上鸟类目标的分布密度、速度和飞行方向。2009年,荷兰和比利时空军提议基于FlySafe-BAM 提供鸟情预警服务,由荷兰皇家气象研究所(Royal Netherlands Meteorological Institute, KNMI)组织完成[32]。当时该服务采用的雷达数据由2 部防空雷达(比利时1 部+荷兰1 部)和5 部气象雷达(比利时3 部+荷兰2 部)提供,能够预测北海、荷兰北部、比利时中部和比利时东部的鸟类迁徙。 随着2013年ENRAM 项目的启动,KNMI 通过与波兰、荷兰、比利时等国空军及德国联邦国防军的合作,雷达数据源进一步丰富,FlySafe-BAM 的功能得到进一步提升,能够对荷兰、比利时和德国境内的鸟情进行准实时预报,以鸟情信息通告(Bird notice to AirMen, BirdTAM)的形式发布给飞行员[33]。目前,FlySafe-BAM 中的鸟情检测算法仍在持续改进,随着加入飞行安全计划的成员国数量的增加,FlySafe-BAM 信息服务的范围将继续扩大。

2 气象雷达飞鸟目标回波特性

本节从目标回波幅度、目标高度分布、飞行速度和方向等方面分析气象雷达获取的飞鸟目标回波特性,并讨论其主要特点及与专业探鸟雷达的区别,为基于气象雷达的鸟情信息提取奠定基础。

2.1 目标回波幅度

比较降水与迁徙鸟群的雷达回波幅度具有一定的指导意义。气象学上所谓的Z-R关系反映了雷达反射率Z与降雨强度R之间的关系,通常表示为[34]

式中:Z的单位为mm6/m3,R的单位为mm/h,且ZdB=10lgZ。

反射率与空气中液态水含量之间的关系可由经验公式表示为[35]

式中:M表示单位体积内的液态水含量, kg/km3。在中等降雨情况下(R=1 mm/h,Z=200 mm6/m3)空气液态水含量为M=7×104kg/km3。迁徙鸟群在迁徙高峰期的分布密度可达到100 只/km3,按照雀形目单只鸟相当于50 ml 的水球计算,其液态水的含量能达到5 kg/km3,远低于中等降水[36]。

雷达后向散射不仅取决于散射体的质量和密度,而且与散射体的形状密切相关。瑞利散射条件下,降水密度为5 kg/km3所对应的反射率很低(-27 dBZ)。但是,由于鸟类目标的大小与雷达波长相当,其散射特性不服从瑞利散射,而更接近共振散射,因此在C 波段,100 只鸟/km3会产生约3 dBZ 的反射率,但仍远低于中等降水的反射率。图3 给出了2 次实验获取的鸟类目标分布密度与雷达反射率的对应关系,其中实线为鸟类目标分布密度的均值,带颜色的区域为标准方差分布[37]。2 次实验均采用了OPERA 系统中的C波段气象雷达,2 次实验由于实验地点和时间的不同,存在一定的季节性差异,气象雷达探测到的是不同的鸟群目标,导致相同雷达反射率时,实验2 对应的鸟类目标分布密度略低。

图3 鸟类目标分布密度与雷达反射率的对应关系[37]Fig.3 Relationship between bird target distribution density and radar reflectivity[37]

2.2 目标高度分布

鸟类迁徙高度一般低于1 000 m,小型鸣禽的迁徙高度不超过300 m,大型鸟可达3 000~6 300 m,个别种类可以飞越9 000 m,鸟类夜间迁徙的高度往往低于白天。候鸟迁徙的高度亦与天气有关[38]。天晴时,鸟飞行较高;在有云雾或强劲的逆风时,则降至低空飞行。

欧洲的研究者利用气象雷达,对秋季的迁徙鸟群进行了观测,收集的数据包含了鸟群在多个夜晚的大规模迁徙活动,并采用专业探鸟雷达做参照验证。为避免地杂波的干扰,观测只针对200 m 以上的空域,气象雷达和专业探鸟雷达均能实现较好的覆盖。该研究使用迁移流量率(Migration Traffic Rate, MTR)来比较不同雷达系统观测到的迁移强度,MTR 表示1 h 内通过垂直于迁徙方向1 km 的虚拟样带的鸟类数量,是一种结合鸟类分布密度和飞行速度的流量测量方法,能够反映通过给定区域的鸟类数量[39]。图4比较了3 种雷达系统获取的不同高度上的鸟类MTR 分布[40]。气象雷达和探鸟雷达在9月5日—10月31日整个迁徙季的数据进行了对比,并增加了10月5日—10月16日之间的垂直扫描海事雷达的数据。总体上看,这3 类雷达系统获取的MTR 高度分布数据具有较高的相关性。但是,在较低的高度上(200~400 m),专业探鸟雷达探测到的鸟类目标数量要明显高于气象雷达和垂直扫描的海事雷达。可见,气象雷达在低空探测中易受地面杂波干扰,其更适于探测具有一定飞行高度的迁徙鸟群。

图4 鸟类迁徙中的飞行高度分布[40]Fig.4 Flight height distribution in bird migration[40]

2.3 目标飞行速度和方向

迁徙鸟速度大多在30~70 km/h,迁徙中每天飞行6~8 h。研究表明,候鸟的迁徙速度受气流的影响,顺风快、逆风慢,同时也受气温和季节的影响[41],故而不少鸟类迁徙多选择在上升气流活动强烈的白天或季风时节。这点在猛禽迁徙中表现尤为明显,其经常成群结队以盘旋滑翔方式向前方作滚动式迁徙。

一般情况下,气象雷达测量的鸟类飞行速度要略低于专业探鸟雷达,原因在于飞行速度的计算是基于测量单元内包含的所有鸟类的径向速度,导致气象雷达估计的飞行速度可能会低于真实值[42]。只有单元内所有飞鸟目标的飞行方向保持一致,才能获取更接近真实值的平均速度估计值。在较低的高度上会发生更多的非迁移运动,并且鸟类的运动会受到地形的影响,这也会导致气象雷达在低空观测到的目标速度较低。同时,在低海拔地区,少量杂波也可能与相对较弱的鸟类目标回波混合,这可以解释气象雷达探测到的一些较低的目标速度。可以认为,当定向散射很小时,由气象雷达得到的低空平均地面速度是可靠的。图5 为基于OPERA 系统中的C 波段气象雷达获取的某日飞鸟目标飞行速度与飞行方向的分布统计结果,并与探鸟雷达的探测结果进行对比验证[37]。在鸟类目标飞行方向的测量方面,气象雷达与专业探鸟雷达的探测结果在相互验证的同时也存在一定差异。尤其是在低空,气象雷达探测的鸟类飞行方向比专业探鸟雷达的探测结果更为集中,这是由于气象雷达给出的目标飞行方向是基于每个扫描单元中所有目标的运动方向均值得出的,而专业探鸟雷达能够给出每个目标的运动方向。当然,在鸟类大规模迁徙的夜晚,目标的运动方向更趋一致,气象雷达与专业探鸟雷达的探测结果也更为吻合。

图5 鸟类目标的飞行速度与方向分布示例[37]Fig.5 Example of flight speed and direction distribution of bird targets[37]

2.4 气象雷达探鸟数据的主要特点

一般情况下,气象雷达开展鸟情探测需要相对较长的观测时间以获取足够的观测数据,进而提取出鸟类迁徙的时间段、迁徙目标空间分布、迁徙方向、飞行高度分布以及平均速度等信息。在迁徙时间方面,由于大部分气象雷达的全空间扫描周期较长(>3 min),存在迁徙飞鸟错过雷达空间扫描空域的可能性,因此需要采用基于鸟情大数据的统计手段获取飞鸟的迁徙时间信息,通常采用时段进行描述,时间分辨率介于30 min~1 h之间,其精度与观测样本数据量正相关;在迁徙目标空间分布方面,气象雷达的距离和方位分辨率通常较低(距离分辨率>100 m,角度分辨率≈3°),因此,基于气象雷达标定的鸟群位置信息通常取决于雷达体扫描单元,且距离越远,雷达体扫描单元越大,水平分辨率在几百到上千平方米之间;在迁徙方向和飞行速度方面,由于雷达只能测量目标的径向速度,且鸟类迁徙过程中飞行方向存在一定的不确定性,因此其飞行方向及速度通常是基于气象模型和多方位观测数据进行反演获得的,分辨率较低且存在一定误差,仅能实现量级上的估计。表3 归纳了气象雷达探鸟数据的主要特点。

表3 气象雷达探鸟数据的主要特点Table 3 Main characteristics of weather radar for bird data

3 气象雷达鸟情信息提取技术

基于气象雷达数据提取鸟情信息,需要首先排除杂波和气象信息的干扰,进而针对飞鸟目标的回波特征提出一系列鉴别指标。本文在总结当前研究成果的基础上,介绍了基于鸟情大数据和深度学习的鸟类迁徙规律分析方法,并对当前气象雷达鸟情信息提取技术的难点和相关解决方法进行了梳理。

3.1 杂波抑制

杂波是多普勒雷达进行反射率和平均速度观测时需要剔除的信号。在低空空域,杂波可能是由地面的后向散射或雷达周围高层建筑的局部波束遮挡引起的[43]。由于地面产生的天线旁瓣功率散射,通常会检测到沿方位方向的杂波线。典型的地杂波在近距离处最强,并且随着海拔高度的增加而显著减小。大气中的逆温会导致雷达波束的异常传播,使波束向下折射到地球表面。这可能会导致在更大范围内从地面产生强散射。鸟类密度估计是基于晴空观测条件下的鸟类目标回波信号,这种回波的强度往往较弱。因此,鸟类反射率估计对杂波引起的伪回波污染非常敏感。除目标反射率的观测之外,杂波还会影响目标平均径向速度观测的质量,这在鸟类的检测中非常重要。相关研究表明,叠加在降水气象信号上的地杂波会导致信号质量降低,并使平均径向速度向零偏移[44]。本节将对已经开发的气象雷达鸟类目标检测算法中使用的杂波抑制策略进行讨论。

多普勒信号中的地杂波干扰可以通过时域或频域的数字滤波进行抑制。在多普勒频谱中,地杂波在零频率附近产生一个窄的峰值。这种低频成分可以通过应用一个尖锐的高通滤波器从多普勒信号中去除,从而产生剔除杂波的径向速度场和反射率场。在鸟类迁徙过程中,多普勒杂波滤波的效果不如晴空条件。通常情况下,地杂波强度约为0~15 dBZ,而鸟的回波约为-2~5 dBZ[45]。接近雷达的杂波回波与鸟回波具有相似或更大的幅度,这意味着杂波比信号的强度至少高出几分贝。在鸟类迁徙过程中,动态多普勒杂波滤波的效果明显降低。

在非鸟类迁徙季节的几个晴天内,通过对所有距离门进行平均反射率计算生成静态杂波图,是一种普遍采用的杂波抑制方法[46]。所有反射率平均值高于-10 dBZ 的距离门将被永久剔除。由于地杂波在多普勒频谱的零频率附近产生窄的峰值,可以排除多普勒速度在-1~1 m/s 间隔内的所有距离门。海杂波往往具有非零多普勒速度,因此,临海的气象雷达的杂波抑制难度会更大。

3.2 气象信息剔除

在鸟类探测的过程中,降水是另一个需要充分考虑的“杂波”来源。与鸟类目标的散射相比,降水的回波强度通常更大[47]。被误认为鸟类的少量降水信号可能会导致对鸟类目标密度的高估,因此剔除所有非鸟类目标回波的区域至关重要。KNMI 的气象学家通常使用7 dBZ 的低反射率阈值来区分降水信号和低反射率的晴空回波[39]。在C 波段,典型的鸟类反射率一般不足7 dBZ,而各类气象学上不显著的降水信号也是如此,需要利用多普勒径向速度进行区分[48]。

针对每个仰角的扫描结果,通过搜索每个单元筛查出反射率高于某阈值的区域。该阈值需要设置的足够低,以确保包含所有降水区域,但此时不可避免地会选择部分鸟类迁徙区域,需要附加分析步骤对降水单元加以确认。设置反射率阈值为0 dBZ,能够包含观测到的大多数降水气象信号。虽然在许多情况下,当选择低于0 dBZ 的阈值时,降水滤波更有效,但在部分情况下,其可能导致降水和鸟类迁徙区域合并为一个单元,从而降低气象剔除的准确性。因此,对径向速度的整体分布特征分析对于随后消除这些降水回波非常必要[49]。降水区域的单元格通常是均匀填充,但鸟类目标的回波区域一般不满足均匀填充。对于筛选出的降水单元,要求其与至少5 个高于反射率阈值的其他单元相邻,这保证了只选择出均匀填充的单元格。

反射单元的径向速度特性也可以作为区分降水和迁徙鸟类目标回波的标准之一。从径向速度场导出的纹理场等于局部径向速度的标准差,而每个距离门都需要考虑所有直接相邻的单元。径向速度Vr的局部标准偏差σv为[37]

图6 绘制了气象雷达分别在较强的鸟类迁徙和对流降水期间探测的单元平均局部标准偏差σv和反射率Z的对应关系,其中红色点为迁徙鸟目标回波,蓝色点为降水回波[37]。图中所示数据由2 个实验地点分别部署的2 部气象雷达在秋季迁徙季获取。通过实验现场验证,鸟类迁徙过程中未出现降水。由图6 数据可知,迁徙鸟回波数据满足σv>5 m/s 和Z<15 dBZ 这2 个条件。显然,图6(a)中的回波数据较图6(b)更易区分,这是由于前者采用的气象雷达的距离分辨率更高。由于实验1 中采样单元的体积较小,获取的数据精确度更高,计算出的反射率结果也更为准确。实验2 中采用的气象雷达采样单元相对较大,也会导致径向速度在更多鸟类目标上进行平均,从而降低径向速度标准差中的纹理值。此外,实验2观测区域内较强的地杂波可能导致径向速度的去混叠误差,进而增加降水回波的径向速度标准差。

图6 迁徙鸟群与降水信号的特征分布[37]Fig.6 Characteristic distribution of migratory birds and precipitation signals[37]

3.3 飞鸟目标特征提取

对于单极化气象雷达来说,3.2 节所述的目标反射率Z和局部标准偏差σv是识别飞鸟目标的主要特征依据。其中,局部标准偏差σv也称为谱宽度,其是一个体扫描单元中所有散射子径向速度的统计标准差,描述了体扫描单元内全体散射子运动速度的一致性。随着双极化气象雷达的日益普及,其能够提供更为丰富的目标回波数据,包括差分反射率ZDR、差分相位ψDP和极化相关系数ρHV等,为飞鸟目标特征提取与识别创造了更好的条件[50]。

在双极化条件下,反射率既可以测量水平波的贡献,也可以测量垂直波的贡献,并且通常与每个极化方向上散射体的大小有关。利用对数变换规则,差异反射率ZDR可以定义为水平极化ZH和垂直极化ZV下测量的雷达目标反射率之差[51]:

对于气象散射体和波长远大于散射体尺寸的散射体,ZDR值与被测散射体的物理形状(即长宽比)直接相关,因此ZDR的大正值表示水平方向的扁圆物体,大负值表示垂直方向的长形物体,接近零的值表示球体或随机取向粒子的集合。但是,这种简单的解释通常不适用于生物散射体,存在一定偏差。对于鸟类这样的非球形散射体,每个极化方向都会受到谐振影响。对处于谐振区的飞鸟目标,稍微增大鸟的尺寸可能导致垂直极化RCS 的增加和水平极化RCS 的减少,从而产生较低的ZDR值。显然,对于飞鸟目标来说,ZDR和目标长宽比之间并不存在简单的对应关系。目标RCS 不仅取决于其尺寸,还与雷达波长有关。为避免相邻雷达站点之间的无线电干扰,NEXRAD 网络中相邻站点的雷达波长设置也略有差异。Melnikov 等的研究表明,波长的轻微变化(10.0 cm 和11.1 cm)会导致同一个鸟类散射体的ZDR值从+10~-5 dB 的变化[52]。

差分相位ψDP是水平和垂直极化之间的总测量相位差,包括发射时的初始差分相位ψt、电磁波发射过程中的相位偏移0.5ΦDP、目标散射造成的相位偏移δ、电磁波返回过程中的相位偏移0.5ΦDP以及接收回波时的相位偏移ψr,表示为[53]

式(5)的结果是一个介于0°~360°之间的值。在气象应用中,该值被解释为水平波滞后于垂直波的程度。在生物目标探测中,垂直波可能滞后于水平波,导致负相位偏移。实际上,在晴空条件下探测生物目标,传输过程中造成的相位偏移(ΦDP)几乎可以忽略,总的差分相位可以减少到雷达系统的贡献(ψt和ψr的组合)和散射体引起的相移(δ)。为充分考虑雷达系统的贡献,大部分雷达系统都会进行系统相位校准,例如,NEXRAD 雷达系统在应用中会设置一个60°的差分相位作为所有站点的校正标准。值得注意的是,许多实时雷达可视化工具通常不能给出差分相位(°),而是给出距离导数((°)/km),原因在于此类产品是专门为气象应用设计的,因此不能针对特定生物目标提供任何特征数据。

雷达接收机记录的回波是大气传播、生物目标散射和雷达系统硬件本身变化的复杂组合,可能导致在每种极化条件下被测信号的振幅和相位受到不同的影响,因此,可以通过计算二者之间的相关性进行比较,通常称为极化相关系数ρHV。对于每个分辨单元,这种相关系数的大小取决于2 种极化跨越多个脉冲的接收信号的相似性。对于含有固定散射体的采样单元,不存在由散射体引起的脉冲间变化,ρHV=1。如果散射体的运动、形状或方向更加多变,则每种极化信号的相关性降低,ρHV值将变小。与雨滴不同的是,鸟群等生物目标的形状、身体位置、方向和空间分布都是高度可变的,直接导致ρHV值降低。此外,鸟类飞行中通过翅膀拍打和弯曲导致机体几何结构的实质性变化,也会改变其后向散射特性[54]。当然,这些散射机制及其随动物体位的变化规律仍待研究。

基于以上气象雷达目标特征数据,能够生成空间特征纹理图像。特定雷达产品的某个像素位置的纹理通常由其周围像素的某种可变性统计值来定义的。此类计算中包括的像素可能在距离和方位角上跨越不同区域。极化纹理与生物目标在分辨单元内及其之间的特性和行为的可变性有关。例如,当采样单元的组成不均匀时,相邻单元将具有不同的极化特征,导致较高的空间标准差。如果目标行为方式不同,在同类目标组成的情况下仍然会出现较高的标准差。例如,2 个相邻单元可能都包含相同的鸟类物种,但如果观测角度不同,仍会产生不同的ZDR值和较高的ZDR纹理。图7 为NEXRAD 网络中位于田纳西州的一部气象雷达于2015年5月3日获取的夜间迁徙鸟群目标回波并生成的多个指标数据产品,包括反射率Z、径向速度Vr、谱宽度σv、差分发射率ZDR、差分相位ψDP以及相关系数ρHV,该部雷达此时的仰角为0.5°,图中所示的探测半径为300 km。显然,ZDR随方位变化显著,负值向西呈楔形延伸,正值向东南呈楔形延伸[8]。实际上,以上区域内不太可能存在鸟群飞行行为的空间变化,ZDR 的变化纯粹是由雷达不同的观测角度造成的。此外,ZDR 的这些负值也并不表示鸟类目标的垂直尺寸更大。基于多个特征指标判断,该迁徙鸟群正在朝东北方向运动。

图7 鸟群目标特征指标纹理图[8]Fig.7 Texture map of bird flock target characteristic index[8]

实际上,鸟类在迁徙过程中的飞行方向和速度存在一定的不确定性,且不同种群的鸟群在迁徙模式上存在差异,导致气象雷达对鸟类迁徙方向及飞行速度的估计精度不高。此外,由于雷达仅能根据回波的多普勒特征获取目标的径向速度信息,实际飞行方向和速度未知,因此也增加了对迁徙鸟群目标的飞行方向和速度进行估计的难度。目前主要的解决方法包括:研究基于气象雷达数据的风向与风速评估模型,结合迁徙鸟群的飞行特征以及不同时空窗口上的目标回波信息,对鸟群目标的迁徙方向和速度进行评估,并结合其他观测手段和信息来源加以修正,构建以气象雷达观测数据为核心的鸟群迁徙模式综合评估体系。当然,以上方法获取的目标信息仍然是对迁徙鸟群目标态势的整体把握,无法获得单一目标的精细化特征。

3.4 机器学习技术

每年都有数以十亿计的动物在季节性迁徙期间穿越地球,但由于它们的活动不可预测,导致对其进行预测的难度较大。利用气象雷达提取并积累的海量鸟情信息,通过机器学习技术能够发现其中的规律[55]。牛津大学和康奈尔大学的科学家利用NEXRAD 气象雷达网获取的鸟情信息,基于23年的春季观测资料,采用梯度提升回归树建立了一个大陆尺度的鸟类迁徙预测模型,以确定气温、气压、风力等气象条件与鸟类迁徙强度之间的关联,相关成果发表在2018年的《Science》杂 志 上[56]。该 模 型 解 释 了 美 国0~3 000 m 高度范围内鸟类迁徙强度的81%的变化情况,并且能够提前1~7 天较好地预测迁徙变化(解释了62%~76%的变化)。该模型还发现,在迁徙高峰期,美国各地迁徙过境的鸟类每晚可能超过5 亿只。图8 给出了2008—2017年迁徙高峰期从美国大陆迁徙过境的鸟类数量变化情况,迁徙规模在5月份达到峰值。图9 给出了不同日期和温度条件下迁徙密度的预测值,图9(a)所示热力图中的每个点代表了一部气象雷达在夜间获取的数据,图9(b)为春季的3 个特定日期中鸟类迁徙密度随温度的变化情况。

图8 美国大陆的夜间鸟类迁徙数量[56]Fig.8 Number of nocturnal bird migration in the continental United States[56]

图9 迁徙密度与时间和温度的关系[56]Fig.9 Relationship between migration density and time and temperature[56]

马萨诸塞大学阿默斯特分校的研究者同样利用NEXRAD 气象雷达网获取的鸟情信息,将隐变量模型与期望最大化算法相结合,采用深度学习技术建立更为精准的目标检测器,发现了以前未知的鸟类栖息地,为生物学家了解群居鸟类的活动规律提供了重要信息[57]。研究发现,大量群居鸟类通常在日出前15~30 min 离开栖息地,难以被人类发现。图10(a)第1 行为气象雷达回波图像,方框中的位置为采用以上机器学习技术发现的具有“环形”特征的燕群栖息地,中间一行反映了燕群飞离栖息地过程中的反射率变化情况,最下面一行显示了燕群从栖息地中心飞离的径向速度变化情况,靠近雷达站的飞鸟径向速度为负(绿色),远离雷达站的飞鸟径向速度为正(红色);图10(b)为实地拍摄的燕群图像。

图10 气象雷达获取的鸟类栖息地雷达回波与现场验证图像[57]Fig.10 Radar echo and field verification image of bird habitat obtained by weather radar[57]

目前,现有基于气象雷达的鸟群观测方法基本上仅能够实现飞鸟和其他气象目标以及杂波的区分,且精度仍然存在提升空间,尚不具备区分不同类型鸟群目标的能力。因此,现有气象雷达的鸟群检测方法仅能够初步解决“有无”问题,尚无法精确判定迁徙鸟群的类型并实现对迁徙鸟群的精确描述。目前主要的解决方法包括:从理论建模的角度深入研究不同类型鸟种、密度以及飞行队形对于气象雷达视角下鸟群目标特性的影响,结合探鸟雷达、光电、卫星追踪等其他观测手段和信息来源构建迁徙鸟情目标信息特征库,实现迁徙鸟群信息的全方位、多尺度描述。此外,以深度学习为代表的人工智能新技术也可能成为迁徙鸟群目标分类的辅助手段。

3.5 交叉验证技术

基于气象雷达数据的鸟情信息提取算法的正确性,需要人工观测、卫星追踪、专业探鸟雷达等其他鸟情观测技术加以验证[58]。瑞士鸟类研究所利用BS-MR1 型探鸟雷达、垂直扫描的海事雷达以及目标跟踪雷达开展了大量的交叉验证实验,各类探鸟雷达如图11所示[40]。2.2 节 和2.3 节对比了以上探鸟雷达和气象雷达获取的鸟类目标飞行高度、速度和方向等数据,二者在数据分布的统计特征上趋于一致,充分验证了基于气象雷达的鸟情信息提取算法的有效性。

图11(a)所示的BS-MR1 型探鸟雷达由瑞士鸟类研究所研发并专门用于鸟情观测,采用垂直对空凝视的天线体制。该雷达为25 kW 脉冲X波段,并采用长短脉冲相结合的模式,其不仅能获取鸟类目标的速度和飞行方向,还能根据翅膀拍打模式对鸟类进行初步分类。同时,该型雷达能够在考虑不同大小目标的检测概率的情况下,计算出各类目标的MTR 值。在这项研究中,BSMR1 雷达使用4 种工作模式探测回波,分布为静态短脉冲、旋转短脉冲、静态长脉冲和旋转长脉冲。MTR 的计算限制为200~800 m 高度范围内检测到的短脉冲回波,因为短脉冲条件下小鸟的最大检测范围不超过800 m。在800~1 400 m 高度范围内,仅使用长脉冲检测回波。只有在旋转模式下才能检索到鸟类的飞行行为数据。

图11(b)为25 kW 的X 波段垂直扫描雷达,配备2.17 m 的T 形天线,波束宽度为22°,方位角1°,每分钟旋转34 圈。天线的旋转方向与预期的鸟类主要飞行方向保持一致,以便探测鸟类的较长轨迹。由于只进行垂直旋转,鸟在雷达波束内的飞行方向无法确定,限制了某些信息的获取。在数据采集过程中,雷达选择长脉冲的工作模式。为了排除昆虫目标,数据处理中忽略了短于200 m且连续回波<4 次的轨迹,还排除了距雷达300 m以内的轨迹以及速度<30 km/h或>100 km/h的轨迹。MTR 和鸟类飞行速度的计算都基于所有鸟类都穿过平行于雷达旋转轴的波束的假设。因此,对于飞行方向偏离雷达旋转轴的鸟类,可能低估其轨迹长度和飞行速度。

图11 用于交叉验证的各类探鸟雷达[40]Fig.11 Various avian radars for cross validation[40]

图11(c)所示的跟踪雷达为200 kW 的X 波段雷达,波束宽度为1.5°,通过手工操作进行目标跟踪定位,每次只能跟踪一个飞鸟目标。在跟踪目标的精确位置时,每秒记录一次数据,对目标的飞行高度、速度和航迹方向进行精确测量。在单个目标情况下,根据翅膀扇动模式的回波特征,识别出非鸟目标、鸟目标、雀形目或鸟群。该跟踪雷达被用于目标轨迹方向和速度的交叉验证。由于手动选择跟踪目标以及目标跟踪的持续时间可能导致跟踪目标数量的偏差,因此,跟踪雷达数据不用于估计迁移强度和海拔分布的验证。需要强调的是,气象雷达是针对气象类的大尺度目标进行全空间观测,因而通常时空分辨率较低,获取的目标信息量较少。因此,与专业探鸟雷达相比,气象雷达的扫描速率较低,空间分辨率不足。通常,大部分迁徙鸟群目标仅分布在气象雷达的一个体扫描分辨单元内且时空不确定性均较高,导致气象雷达在全空间扫描过程中可能无法扫描到迁徙鸟群,获取的有效信息不足。目前主要的解决方法包括:延长观测时间以获取足够数量的观测样本,通过大数据方法实现对迁徙鸟情态势的统计学特征描述。

4 基于气象雷达的鸟类学研究与应用

基于气象雷达组网的鸟类探测技术在大陆尺度鸟情监视预警方面具有独特的优势,日益成为鸟类生态学研究与军民航鸟击防范的重要技术支撑。

4.1 鸟类生态学

长期以来,由于缺乏大尺度鸟情观测的技术手段,掌握夜间鸟类迁徙的飞行路线一直是一个科学难题。欧洲的研究者利用OPERA 气象雷达网络调查了鸟类在欧洲范围内夜间的飞行路线,这是首次利用欧洲天气雷达网络研究大陆尺度的大规模鸟类迁徙活动[59]。图12 为2016年秋季迁徙季节,基于从斯堪的纳维亚半岛北部到葡萄牙的70 个气象雷达站提取的生物信息,绘制的鸟类迁徙主体方向图,显示出欧洲部分地区的迁徙强度[60]。图12 中每个监测位置的平均迁移方向用黑色箭头表示,每个位置的平均MTR 用圆圈的大小和颜色表示,取整个采样期的平均值(2016年9月19日—10月8日,仅包括夜间数据)。每个箭头的宽度表示迁徙方向的集中性,较宽的箭头表示更集中的迁移方向,较细的箭头表示更分散的方向分布。平均而言,基于所有监测位置的20 个夜间采集的雷达数据,每千米样带每小时有389 只鸟通过。在迁徙强度最高的夜晚,每个雷达站监测到的MTR 平均值为1 621 只/(km·h),但迁徙强度在地理位置和时间分布上有很大差异。总体的迁移方向为西南向,法国中部的强度最高。研究发现,迁移行为与风力条件密切相关。通常在顺风条件下,风速每提高1 m/s,MTR 值将相应提高约33只/(km·h。)

图12 欧洲部分地区的季节性鸟类迁徙[60]Fig.12 Seasonal bird migration in parts of Europe[60]

4.2 鸟击防范

AHAS 是基于美国空军鸟击危险减少计划、鸟 类 躲 避 模 型(Bird Avoidance Model, BAM)以 及NEXRAD 气 象 雷 达 网 研 发 的[61]。其中,美国空军鸟击危险减少计划始于20 世纪70年代,最初致力于机场对于鸟禽的主动和被动管理,20 世纪80年代中期以来开始实施低空飞行任务鸟类躲避计划,侧重历史鸟击事件数据的建模分析和实时的鸟类躲避研究。美国空军BAM 建设始于1983年,该模型在1998年实现了地理信息系统与雷达数据的结合,能够提供 全 美 鸟 类 的 分 布 和 数 量 信 息[62]。 美 国NEXRAD 国家气象雷达网,拥有160 部气象雷达,是唯一能够覆盖美国全境的探测网络,每6~10 min 更新一次数据。

AHAS 利用气象雷达网提供的数据,采用神经网络模型等识别算法,实时从气象雷达数据中提取鸟类活动信息;系统利用雷达探测数据,以及对历史鸟情数据的积累分析,对特定区域鸟情危险程度进行评估,风险等级分为低、中、高3 级,评估信息可与地理信息系统结合显示,也可通过邮件或短信的方式向用户发出预警。AHAS 通过2 个因素评估鸟击风险:“严重程度”基于雷达回波强度得出,“鸟击概率”取决于鸟群活动的多边形区域所占百分比,二者相乘得出鸟击风险值。AHAS 的优势是基于现有国家气象雷达网覆盖美国全境,适于探测成群飞行、沿海岸线飞行、或大规模迁徙飞行的鸟类,并能够在地理信息系统中标识,直观显示某地区或某航路的鸟情态势分布,便于空军提前做好鸟防准备或调整训练安排。美国空军自2003年使用AHAS 以来,其鸟击事件数量及损失下降了50%以上,鸟击防范效果显著。图13(a)[61]为AHAS全美鸟情动态预警,其中红色为高度危险、黄色为中度危险、绿色为低度危险、蓝色为气象信号。

图13 美欧鸟情分布态势显示界面Fig.13 Display interface of bird distribution in US and Europe

欧洲基于OPERA 气象雷达网络的数据和鸟击风险评估模型FlySafe-BAM,能够实现对欧洲部分地区的BirdTAM 信息发布[63]。图13(b)[63]为德国、荷兰和比利时范围内,雷达合成的鸟类和昆虫等生物目标在约1 km 以下的低空分布情况,由BirdTAM 值表示鸟类或昆虫的密度。BirdTAM 值遵循0~8 的对数标度,并扩展为≥9 的值,以区分大规模鸟类迁徙和降雨信息。鸟和昆虫的回波显示为绿色、黄色、橙色、红色,有时甚至更暗的颜色;较大的降水范围由黑色或深色区域表示;山脉等地面杂波呈现为小的黑暗区域,如德国南部和东南部以及比利时东部。

5 中国气象雷达鸟情预警系统建设构想

在中国建立覆盖全国的鸟情预警系统有助于实时掌握鸟类迁徙信息,及时向机场、航空公司、空管等相关单位通报高风险鸟情,从而做好协同防范,提高鸟击防范能力。本节将从探测覆盖范围和雷达系统性能2 个方面,讨论在中国建立基于气象雷达组网的全国鸟情预警系统的初步构想。

5.1 探测覆盖范围

中国从20 世纪90年代后期开始建设的新一代天气雷达网,在灾害性天气监测和预警服务方面发挥了重要作用,同时带动了全国气象雷达产业、技术和人才的快速发展,取得了突出的社会经济效益。

2017年6月,中国气象局发布了《气象雷达发展专项规划(2017—2020年)》(以下简称《规划》)。《规划》指出,截至2016年底,全国已经完成233 部新一代天气雷达建设;中国气象局统筹建设的X 波段天气雷达共42 部,由地方自主建设的X 波段天气雷达约200 部;完成部分天气雷达的双偏振升级改造;天气雷达近地面1 km 的覆 盖 范 围 约2.2×106km2[64]。世 界8 条 候 鸟 迁徙路线中,有3 条经过中国:第1 条是西太平洋路线,经过中国东部沿海省份;第2 条是东亚澳洲的迁徙路线,经过中国中部省份;第3 条是中亚、印度的迁徙路线,翻越喜马拉雅山,经过青藏高原等西部地区。其中,东亚澳洲候鸟迁徙路线是全球8 大路线中物种数量最多的区域,在中国,这条候鸟迁徙路线从东北的黑龙江到黄渤海沿岸,南至广东南岭等地,还覆盖内陆湿地鄱阳湖、洞庭湖等[65]。按照相关规划,新一代天气雷达网已经覆盖了中国大部分机场和候鸟主要迁徙路线,能够初步满足全国机场的鸟情预警需求,且这一网络还在不断地完善和升级过程中。目前,按照统一技术标准,已完成112 部新一代天气雷达进行技术升级,并使用X 波段天气雷达对新一代天气雷达网的探测盲区进行补充。

5.2 雷达系统性能

中国的新一代多普勒天气雷达是基于美国WSR-88D 天气雷达先进技术的基础上研发的,具有24 h 无人职守、实时监控、实时标校、高精度、高可靠性等特点,其技术性能已达到世界领先水平,可实现雷达联网、数据互传和资料共享。

近年来,国外气象雷达的整体发展趋势具有从单偏振(极化)探测向双偏振探测发展的特点,美国已于2012—2013年完成了全国天气雷达网的双偏振技术改造并投入业务化应用。双偏振天气雷达与目前常用的单偏振天气雷达相比,能够获取目标的形状、尺寸大小、相态分布、空间取向等更为详细的信息,有助于提高雷达目标探测的准确性、精度和探测数据的质控能力。目前,中国已经启动了天气雷达的双偏振升级改造,开展了双偏振雷达业务应用试验,对双偏振雷达探测能力及业务运行情况进行评估,确定双偏振雷达的业务体制、定标方法和技术标准。目前,东部沿海已有51 部雷达升级为双偏振雷达,其中广东省12 部全部升级为双偏振雷达。从技术体制角度来看,中国气象雷达的双偏振技术升级不仅仅为气象目标观测提供了更加丰富的信息,可进一步提升中国的天气预报以及相关气象领域的研究水平,也为基于气象雷达网的迁徙飞鸟目标识别提供了更加丰富的特征信息以及数据支撑。

6 结论与展望

本文在充分研究美国和欧洲气象雷达鸟情预警系统的基础上,讨论了气象雷达鸟类目标回波特性与鸟情信息提取技术,并介绍了基于气象雷达在鸟类学研究与鸟击防范方面的应用情况,最后论证了在中国建立基于气象雷达组网的全国鸟情预警系统的可行性,共得出以下几点结论:

1)气象雷达由于受到技术体制的限制,在探测能力上与专业探鸟雷达有较大差异。主要表现在气象雷达不能获取单只飞鸟目标的确切飞行轨迹,更适于探测具有一定飞行高度的迁徙鸟群,从而在大陆范围内反映出大尺度的鸟类活动态势,因此成为鸟类生态学研究与鸟类迁徙预警的重要工具。

2)中国的新一代多普勒天气雷达与美国的WSR-88D 雷达性能接近,其已经形成且不断完善的气象雷达网的探测范围基本覆盖了全国大部分机场和中国境内主要的候鸟迁徙路线,能够为建立全国鸟情预警系统提供基础雷达数据。

3)基于当前的新一代多普勒天气雷达,特别是经过双偏振升级改造的天气雷达,采用当前的雷达数据处理算法与深度学习理论从雷达回波信息中分离出气象信息和鸟群目标信息,在技术路径上是可行的。

随着气象雷达信号处理技术与硬件系统的不断发展,未来气象雷达探鸟技术呈现出以下几个趋势:

1)精细化。随着气象雷达软硬件系统的不断发展,未来气象雷达将具备更高的空间分辨率和数据更新速率,为在时空维度上实现对迁徙鸟类型以及分布态势的精细化感知创造了条件。

2)组网化。气象雷达的应用领域决定了其具备组网探测的特性,随着气象雷达网络的日趋完善,未来可进一步研究基于气象雷达网的迁徙鸟探测以及动态可视化描述技术,实现区域和全国范围内的迁徙鸟情动态更新。

3)融合化。目前除气象部门,空管及军事部门均配置了专属气象雷达,这些气象雷达与气象部门的专用雷达在探测范围上存在重合或者互补,通过信息融合可实现更大范围内的鸟情观测。此外,空管雷达、专业探鸟雷达以及其他体制的雷达理论上也可以实现不同空域范围内的迁徙鸟探测,且探测精度更高。因此,在日益成熟的雷达数据共享机制的支持下,通过研究不同体制雷达的鸟情信息融合技术,有助于实现全空域范围内迁徙鸟的多尺度时空观测及态势描述,进一步提升军民航鸟击防范能力。

猜你喜欢
鸟群鸟类气象
善于学习的鸟类
气象树
在你灵魂里沉睡的鸟群
《内蒙古气象》征稿简则
我的湿地鸟类朋友
为什么鸟要成群飞翔?
为什么鸟群飞行时不会彼此冲撞?
鸟类
大国气象
美丽的气象奇观