引言

S波段多普勒天气雷达自投入使用以来，便经常在夏秋两季观测到存在明显昼夜变化的晴空回波。随着研究的不断深入，这种昼伏夜出的晴空回波的产生原因被不断揭示。大气湍流造成的折射指数起伏引起的布拉格衍射(又称湍流散射)和昆虫、鸟类的生物活动是大多数这类回波产生的根本原因(张培昌等，2001)。美国国家大气研究中心(NCAR)的Wilson et al(1994)通过对多波段雷达的对比分析，认为昆虫造成的粒子散射是美国佛罗里达州和科罗拉多州地区上空晴空回波的主要原因；魏鸣等(2007)、黄琴等(2018)基于南京地区晴空回波的演变特征及大气垂直结构，认为大气湍流造成的湍流散射是南京地区上空晴空回波的主要原因，并对其偏振信息的产生原因进行解释。虽然目前对这两种不同原因所导致的晴空回波的比例仍未明确，但是不可否认，天气雷达数据对农业、林业和生态等方面研究有着巨大潜力(Bauer et al，2017)。

美国国家海洋和大气管理局(NOAA)国家强风暴实验室的Zrnic and Ryzhkov(1998)开展的试验表明，偏振参量尤其是差分反射率、差分相位和相关系数，在生物回波识别上起着决定性的作用。然而，我国目前天气雷达双偏振改造处于起步阶段，基于偏振参量的识别方法仍有待研究和应用。李丰等(2012)结合地面和卫星资料对非降水气象回波特性进行统计分析，提出了基于模糊逻辑识别和回波分块的非降水气象回波识别算法；黄兴友等(2018)通过对反射率因子进行统计，利用模糊逻辑识别并消除地物回波，但由于生物体造成的粒子散射和大气湍流造成的湍流散射在反射率因子的强度和分布上都较为相似(Hubbert et al，2018)，所以很难用类似的方法识别区分生物造成的生物回波和大气湍流造成的湍流回波。

本文利用生物的定向行为(吴先福等，2006；Chapman et al，2010)，即利用体型较大的迁飞害虫和鸟类飞行过程中因自主飞行和调整迁飞方向进而导致与背景风速的差异，通过滤波的方式筛选出存在速度差值区域并估算该区域为生物迁飞活动的概率，从而实现了利用天气雷达对定向迁飞生物的实时监测。

1 研究区域与资料选择

我国迁飞性害虫主要的迁徙通道分为湘桂走廊和渤海湾通道：在春季害虫随盛行西南风经湘桂走廊迁飞至江淮和黄淮地区，夏季则经渤海湾向东北地区继续移动(付晓伟，2015)。一方面，相较于陆地，海面平坦，地形因素所导致的风速变化较小；另一方面，中国农业科学院植物保护研究所长岛有害生物防治野外科学观测实验站(以下简称农科院长岛实验站)位于渤海海峡中央的山东省长岛县北隍城岛，可提供人工观测资料进行验证。故此，选择渤海湾通道开展有关研究。

分别选用烟台CINRAD/SA和大连CINRAD/SC的S波段天气雷达资料。两部雷达相距158 km，分别距离北隍城岛107和86 km，基本可覆盖整个渤海海峡区域，为监测渤海湾上空迁飞害虫提供了良好的数据支持。

2 识别方法

由于对于气象云雨回波、地物(海)杂波和超折射回波的识别方法国内外已有较为成熟的研究和识别方法(刘娟，2015)，故此，本文识别方法仅区分存在共同定向现象的生物体所产生的生物回波与湍流所产生的回波。为了避免歧义，文中由定向迁飞生物所产生的回波统称为生物回波，而由湍流所产生的回波简称为湍流回波。

2.1 基本假设

算法基于以下两点假设条件：

(1) 湍流的多普勒速度完全等于真实风速；生物运动速度不完全等于风速。

湍流是空气按一定规律产生的具有相似性质局地涡旋运动的湍团，其整体运动即为所在区域风的运动情况。

而对存在共同定向现象的生物体，在飞行过程中虽然为节省体力会选择顺风飞行，但由于其在迁飞过程中需要通过调整身体姿态和朝向对前进方向进行定向调整，这便会造成其运动方向与所处风场产生差异，故此其运动矢量不完全等于实际的风矢量，例如昆虫主动调整头部指向后与所处风场的风向产生夹角，进而使得方向与速度产生偏差。

(2) 不考虑下垫面的影响，雷达照射体积内存在生物的回波强度必大于同一物理性质湍团内不存在生物的回波强度。

湍流的本质依然是空气，而定向迁飞生物在迁飞过程中本质上其运动于湍流中。所以当生物定向迁飞活动出现时，回波信号实际上等同于湍流信号与定向迁飞生物信号的叠加。故此推论，有效照射体积内存在定向迁飞生物的回波强度必大于有效照射体积内仅存在湍流所产生的回波。

2.2 识别基本原理

由于雷达探测到的风场为径向风场，故此在空间中某一点的径向风速可以写为矢量风在雷达径向上的投影乘以风向与雷达径向夹角的余弦。若风场在一定范围内均匀，则雷达所探测到的定向迁飞生物的径向速度可写为

$ {v_{{\rm{migration\_echo}}}} = {\mathit{\boldsymbol{V}}_{{\rm{wind }}}} \cdot \cos \theta + {\mathit{\boldsymbol{V}}_{{\rm{migration }}}} \cdot \;\cos \varphi \pm \Delta v $

(1)

式中：V_wind为背景风场，V_migration为定向迁飞生物自主运动的净速度，θ为回波所在空间位置上风向与雷达径向夹角，φ为回波所在空间位置上生物定向运动方向与雷达径向夹角，Δv为系统误差，范围为0~1 m·s^-1。

由于湍流的运动速度等同于风速，故此当回波为湍流产生时，其速度可写为

$ {v_{{\rm{unmigration\_echo}}}} = {\mathit{\boldsymbol{V}}_{{\rm{wind }}}} \cdot \cos \theta \pm \Delta v $

(2)

当雷达探测到未知区域的径向速度为v_radar时，若已知存在系统误差的背景风场径向速度V_wind·cosθ±Δv，则将其相减所获得速度差v_d为

$ {v_{\rm{d}}} = {v_{{\rm{radar }}}} - {\mathit{\boldsymbol{V}}_{{\rm{wind }}}} \cdot \cos \theta \pm \Delta v $

(3)

若该未知区域由定向迁飞生物所产生，即v_{migration_echo}=v_radar，则将式(1)代入式(3)中，可得

$ {v_{\rm{d}}} = {\mathit{\boldsymbol{V}}_{{\rm{migration }}}} \cdot \cos \varphi \pm 2\Delta v $

(4)

若该未知区域由湍流所产生，即v_{unmigration_echo}=v_radar，则将式(2)代入式(3)中，可得

$ {v_{\rm{d}}} = \pm 2\Delta v $

(5)

由于系统误差Δv最大为1 m·s^-1，所以速度差v_d数值大于2 m·s^-1时，其必然存在V_migration·cosφ项，进而表明该区域存在定向迁飞生物活动。

但在实际过程中，存在阵风或者湍流扰动等因素，所以风场可能偶然存在个别区域的速度差v_d大于2 m·s^-1。同样的，生物回波区域虽然也存在一些速度差小于2 m·s^-1的区域。根据假设(2)，可依照反射率因子的大小对回波进行分块。设A事件为晴空回波的回波块内为存在定向迁飞生物，事件B为回波块内不存在定向生物迁飞，即回波仅由湍流产生。假定回波块内是否存在定向迁飞生物的可能性相同

$ P({\rm{A}}) = P({\rm{B}}) $

(6)

若阵风对背景风的附加径向速度为v_g，则由于定向迁飞生物造成该区域的速度差v_d>2 m·s^-1的概率为

$ \begin{array}{l} P\left({{v_{\rm{d}}} > 2{\rm{m}} \cdot {{\rm{s}}^{ - 1}}\mid {\rm{A}}} \right) = P\left({{\mathit{\boldsymbol{V}}_{{\rm{migration }}}} \cdot \cos \varphi > } \right.\\ \;\;\;\;\;\;\;\;\;\;\;2{\rm{m}} \cdot {{\rm{s}}^{ - 1}}\left. { \cup {v_{\rm{g}}} > 2{\rm{m}} \cdot {{\rm{s}}^{ - 1}}} \right) \end{array} $

(7)

而湍流所产生的回波造成该区域的速度差v_d大于2 m·s^-1的概率为

$ P\left({{v_{\rm{d}}} > 2{\rm{m}} \cdot {{\rm{s}}^{ - 1}}\mid {\rm{B}}} \right) = P\left({{v_{\rm{g}}} > 2{\rm{m}} \cdot {{\rm{s}}^{ - 1}}} \right) $

(8)

则根据式(6)~式(8)式可得

$ P\left({A\;|{v_{\rm{d}}} > 2{\rm{m}} \cdot {{\rm{s}}^{ - 1}}} \right) = P\left({B\;\mid {v_{\rm{d}}} > 2{\rm{m}} \cdot {{\rm{s}}^{ - 1}}} \right) $

(9)

同时，因定向迁飞生物的速度通常远大于2 m·s^-1(Westbrook，2008)，且阵风的空间尺度较小，则针对某一块回波整体而言，P(V_migration·cosφ>2 m·s^-1)远大于P(v_g>2 m·s^-1)，故当v_d>2 m·s^-1，事件A发生的概率远大于事件B。这就意味着，存在定向迁飞生物的回波区域中出现速度差v_d大于2 m·s^-1的概率明显大于不存在定向迁飞生物的回波区域，且概率越大，其为生物回波的可能性就越高。

需要注意的是，当定向迁飞生物方向与雷达径向夹角φ接近±90°时，将会导致V_migration·cosφ的取值快速减小，这就意味着

$ {\lim _{\varphi \to \pm {\rm{ \mathsf{ π} }}/2}}{\mathit{\boldsymbol{V}}_{{\rm{migration }}}} \cdot \cos \varphi \pm 2\Delta v = \pm 2\Delta v $

(10)

此时，无论未知区域速度信息是否包含定向迁飞生物，v_d恒等于±2Δv，算法将在此时失效。

不过，若多普勒速度最小分辨率为0.1 m·s^-1，则生物定向运动速度为2 m·s^-1时，通过计算可得雷达能够分辨的2 m·s^-1速度差异的最大夹角为87.1°，从理论上而言出现算法失效状况的概率较小。但不得不提的是，在较大夹角的情况下，生物回波区域中出现速度差v_d大于2 m·s^-1的概率将有所减小，此时将引入一定的误差。

不过值得庆幸的是，渤海海峡作为生物定向迁飞的通道，以夏季北迁与秋季南迁为主。而烟台与大连的两部雷达恰好位于渤海海峡的南北两方，分别位于北隍城岛的南偏东23°(烟台雷达站)与北偏东47°(大连雷达站)方位上，这使得生物在渤海海峡上空的定向迁飞方向与雷达径向间夹角的cosφ较大，一定程度上避免了算法失效的问题。

2.3 速度滤波原理及方法

当风场均匀时，在同一距离上的径向速度即为一频率单一的正弦信号(李瑞义等，2018)，一次雷达锥扫所获得速度为一个周期，其频率可以视为1周期每体扫(若一次锥扫内速度有n个周期，为了方便，则可称为n circle/volum，简写为n cv^-1)。

由于天气形势的变化，风场必然产生各种不均匀的变化，这种不均匀可以视为在相对大尺度的均匀风场上叠加相对小尺度的局地性风场。同理，在不同尺度风场叠加后所获得的径向速度也可视为不同频率的风场径向速度的正弦信号的叠加。风场的尺度越大，局地性越弱，其正弦信号的频率就越小；风场的尺度越小，局地性越强，其正弦信号的频率就越大。

当风场的正弦信号频率为4 cv^-1时，径向速度在方位角0°~360°范围内，数值可出现四次正负周期变化。这意味着，在以雷达为中心的平面直角坐标系中，在一个象限中可以同时出现正负速度交替，而通常而言在自然情况下极少会出现如此强的风切变，故此可认为4 cv^-1的频率变化可以表示通常的自然风场的通带最大频率。虽然可以通过提高滤波器的截止频率来提高对风场的拟合能力，但是也会降低对于局地性风的剔除，进而削弱昆虫回波识别的能力。

图 2为多普勒速度低通滤波后的结果。图 2a中红色折线为业务运行的烟台S波段天气雷达于2017年5月26日22:54(北京时)雷达半径75 km处0.5°仰角的多普勒速度，蓝线为4 cv^-1低通滤波后的多普勒速度。将两者相减后获得图 2b，其概率分布如图 2c所示。从图 2c中可见，速度差大体分布在±2 m·s^-1的范围内，进一步分析可以发现，在未经数据人工质量控制的情况下，73%的速度差数值小于2 m·s^-1，与2.1节中的假设的结果基本一致。当风场中存在定向迁飞生物时，由于定向迁飞生物的运动方向和风向的差异，会造成在雷达的径向速度上的污染，即定向迁飞生物的速度信息污染原有风场速度信息。若定向迁飞生物在空中分布范围的方位角宽度为20°，可视为在某一方位Q±10°范围内存在可使原有速度值产生的偏差(衰减或增益)的污染信息。将人为添加污染信息与原始的径向速度的低通滤波后的结果进行比较，随机重复试验50次可得到图 3。

图 3中的偏差量设定为衰减20%，这相当于在径向风速10 m·s^-1时产生2 m·s^-1的速度偏差。添加污染信息后的低通滤波结果与添加前具有极高的正相关性，相关系数为0.999 4，协方差为85.976 4；将两者相减取均值，则50次随机重复试验的均值为-0.045 2，相减取方差的均值为0.149 7，当污染造成的偏差达到40%时，两者也仅为-0.085 1与0.416 8，具有高度的一致性。这说明，当雷达探测到的大尺度风场的径向速度中存在小尺度的定向迁飞生物信息时，虽然会对径向速度构成污染，但对低通滤波后风场的影响相当有限。

2.4 算法流程

算法流程可分为五个部分：(1)获取雷达数据；(2)雷达数据预处理；(3)低通滤波求解风场速度；(4)获得局地性速度；(5)回波分块计算概率。获取雷达基数据后先进行预处理，使雷达反射率因子和多普勒速度数据的分辨率匹配，剔除无效值。之后，对雷达多普勒速度数据v_radar使用理想低通滤波器进行滤波处理，即对于高于截止频率的速度完全截止，而对于低于截止频率的速度完全无失真保留。接着，通过最小均方差准则对低通滤波后的速度进行放大以获得近似的风场径向多普勒速度数据v_wind。最后，将v_wind与v_radar相减获得速度差v_d。由于这种速度是将大尺度的风场剔除，仅保留小尺度风场的结果，故此也可在此称为局地性速度。

为提高判断的准确率，剔除个别的存在异常值的点，可通过圆形区域均值滤波对局地速度进行散焦模糊处理。之后，根据回波反射率因子的数值对回波进行分块并计算每块回波速度差大于2 m·s^-1区域的比例，即可视为该区域回波内存在定向迁飞生物的概率。

图 4为业务运行的烟台S波段天气雷达于2017年5月26日22:54(北京时)所获数据。图 4c中使用通带最大截止频率4 Hz理想低通滤波器，在过滤掉频率高于4 cv^-1的多普勒速度信息后，多普勒速度与图 4b中原始速度依然保持着较高的一致性，图 4d为局地速度，均值为1.4 m·s^-1。由于晴空回波在距离雷达较远时，回波在空间的分布不均匀使在同一距离方位角0°~360°范围内的数据缺失，从而在频谱上难以准确识别各种频率分量，造成滤波失效。同理，当定向迁飞生物在空间中大面积覆盖时，算法也有失效的风险。

不过，这两个问题在很大程度上并不会出现在现实中。第一，通常夜间晴空回波在50 km距离内分布均匀，晴空回波夏秋两季强盛时甚至可均匀分布于100~150 km范围内，且因误差产生的局地速度大值区没有在强度上对应的区域，所以一定程度上会缓解对算法识别效果本身的影响。第二，由于辽东半岛特殊的细长形状，使得生物迁飞集中于山东半岛北部与辽南一线的渤海海峡区域，而非大面积横跨整个渤海海域，这就为多普勒速度提供了充足的大尺度低频信息，大大降低算法失效的可能性。

3 算法验证及效果 3.1 算法验证

依照部分学者对于晴空回波及生物回波的研究观点，美国俄克拉何马州业务运行雷达KICT(37.64°N、97.44°W)在2012年9月30日夜间(世界时)观测到了一次明显的生物活动(图 5)(江源，2013)。由于生物回波偏振量的数值取决于生物的体型，所以鸟类和昆虫在偏振量上具有明显的不同(Stepanian et al, 2016)。区别于昆虫差分反射率较大而差分相位不超过40°的特点，鸟类的差分反射率通常在-1~3 dB，差分相位可超过100°，可以基本确定回波主体为随风运动的昆虫运动产生。同时在回波个别区域存在鸟类运动(Zrnic and Ryzhkov, 1998; Broeke, 2013)。

算法结果为图 6a所示。由于当雷达在0.5°仰角进行探测时，靠近雷达的晴空回波会受到地物的污染(Richardson et al, 2017)。故此，距离雷达10 km范围内的识别结果被剔除。

将图 6a与6b、6c比较，低通滤波的识别方法与偏振量识别方法结果在空间分布上一致。以差分反射率在-1~6 dB、相关系数在0.3~0.7为鸟类的偏振量分布区间(江源，2013)，在未使用偏振量仅依靠多普勒速度的情况下，算法成功地从作为背景的昆虫回波中筛选出了鸟类活动，准确率为94.01%，造成误差的主要原因在于回波块边缘和离散点状信号的影响。在图 6a中，鸟类活动区域的计算概率普遍高于0.5，而昆虫活动区域的计算概率基本不足0.2，具有较高的区分度。

3.2 重大迁飞事件

2017年5月30日夜间至31日凌晨(北京时)，农科院长岛实验站观测到大量粘虫迁飞状况，根据人工观察，31日凌晨为高峰时段。29日夜间至30日白天受黄淮气旋东移北上影响，渤海海峡及辽东半岛产生降水，30日入夜后转为槽后，大气较为稳定。

30日23：42(图 7a)，烟台雷达上发现在较为明显的高概率回波区域(38°~38.25°N、121°~121.5°E)。6分钟后，即23：48(图 7b)，该区域回波的生物回波概率显著提高，基本可确认其回波为生物定向移动造成。回波整体向北移动，符合生物活动规律。31日01：03(图 7c)，生物回波向北移动至(38.25°~38.5°N、121°~121.5°E)区域。在01：00大连雷达上可见(图 7d)，在同一区域也具有较高的概率。随后，回波继续向北移动，于02：26开始聚集在大连市旅顺口区老铁山国家级自然保护区(图 7e)，直至03：58(图 7f)，回波基本消失，可认为昆虫降落完毕。此次迁飞过程中，昆虫回波的强度约为14 dBz。由于生物体散射特性较为复杂，对于生物的定量化计算方法的研究正在同步开展。

3.3 讨论

在3.2节中，两部雷达同时观测到一次完整的生物(粘虫)迁飞过程。从30日23：42起至31日03：58结束，迁飞全程约4 h。此次迁飞活动的海上运动路径清晰，与观测人员描述的时间相符。飞行高度约为1.6 km附近，符合生物定向迁飞的高度特征(Drake and Farrow, 1985; 邓望喜, 1981)。

需要注意的是，由于在陆地，特别是山区，其对流层的风场尺度较小，极易产生小尺度的自然风，干扰算法识别结果。故此，该算法在陆地上的使用效果，相较于洋面而言有着较大的局限性。这也是目前为何并为给出虫源地信息的主要原因。同时，由于大气运动的复杂性，雷达探测范围内风场不均匀时，算法的识别能力也受到一定的限制。

在进一步的分析后发现存在一些问题。首先，当生物处于多普勒速度为零的区域附近且生物移动方向近似垂直于雷达径向方向，或者生物运动速度较小时，会使算法失效。这就是导致了为何此次定向迁飞过程中，没有获得明确的起飞地的位置。其次，由于本算法未考虑海杂波的干扰，所以在识别的过程中会发现一些无移动过程的偶发的高概率区域。最后，已有研究表明，鸟类回波与昆虫回波的强度较为接近(Zrnic and Ryzhkov, 1998)，不排除在此次粘虫迁飞过程中存在鸟类混杂在其中。故此，需要在未来的研究中更加深入的挖掘、研究湍流回波、鸟类回波和昆虫回波的相关信息，以期更有效地进行区分。

同时发现一些暂无明确解释的现象。第一，在北隍城岛上空，在相当长的时间段内(31日00—02时)出现高生物回波概率区域，这可能是岛屿产生的地物杂波影响，也可能是部分昆虫选择降落岛上造成，这种现象也存在于老铁山保护区海猫岛(鸟岛)上空(31日02—04时)。第二，如图 8所示，在生物回波移动的过程中，发现生物回波分裂为北支和南支两部分，北支部分继续向北移动，而分裂出的南支回波逐渐消失于海上。造成这种原因的生物学机制尚未有明确的解释。虽然有渔民表示第二天海面确有漂浮着的昆虫尸骸，但客观来说，也不能完全排除这是算法存在的某种未知缺陷。

4 结论

本文基于两条基本假设：(1)湍流的多普勒速度完全等于风速，生物运动速度不完全等于风速；(2)不考虑下垫面的影响，当湍流及大气在水平空间的性质基本一致，生物回波必大于湍流强度，进而推论出定向迁飞生物区别于湍流回波的多普勒速度特征，并通过低通滤波的方式成功分离出生物回波区域。

本算法成功分离了美国KICT业务运行雷达在2012年9月30日夜间观测到的生物活动中定向迁徙的鸟类活动。对于2017年5月30日夜间至31日凌晨在农科院长岛实验站观测到大量粘虫迁飞过程，算法成功地从雷达数据中识别出了具体的迁飞过程。

算法也存在着改进的空间，正如2.4节所说，目前算法对雷达观测的天气背景有一定的要求，需要出现较大范围的晴空回波，才可能利用滤波技术得到背景平均风向风速。同时，由于算法要求背景风场需在一定尺度上均匀，所以在应用算法时雷达探测范围内的大气需要较为稳定，没有对流活动，否则算法会将对流云(泡)和降水识别为迁飞生物。而对于迁飞生物，若其运动能力有限运动净速度较小，或其自身运动方向与雷达径向方向近似垂直时，则会明显降低算法的识别能力，造成算法失效。

不能忽视的是，由于目前的客观条件限制，一方面在此次迁飞中存在大量的无法解释的问题，有待更多的生物活动观测结果的验证和生物学的论证；另一方面，对于定向迁飞生物种类、定向生物数量的识别计算工作还有待进一步开展，对于算法的改进有待更多数据支持和进一步研究。值得期待的是，在未来的几年内，随着我国偏振雷达升级改造完成，利用偏振参量对生物迁飞活动的监测能力将有明显提高，以模糊逻辑和机器学习为代表的算法将发挥重大的作用(张林和杨洪平，2018)，这将为我国生物监测提供更多更好的数据。