机载雷达风切变识别算法研究及在机场预报中的应用

魏鸣，张明旭，2，张培昌，郭巍，周生辉，2

(1．南京信息工程大学气象灾害预报预警与评估协同创新中心，江苏南京210044;2．南京信息工程大学中国气象局气溶胶与云降水开放重点实验室，江苏南京210044;3．上海市卫星遥感与测量应用中心，上海201199)

0 引言

随着我国航空事业的发展，减少危险天气对飞行的威胁日益迫切。在愈加频繁的气象危害中，大气的风切变特别是低空风切变对飞行的影响最大，造成的事故所占比例高(黄仪方，2002)。风切变形成机制复杂，持续时间短且强度较大(武虎子等，2013)。雷暴等强对流天气在一定范围内可产生以垂直风为主要特征的综合风切变区域，尤其以雷暴云体的强下沉气流区和冷性外流气流前缘引起的阵风锋更为严重，对飞行危害最大(张培昌等，2001;王丛梅等，2011)。锋面天气可产生以水平风的水平切变和垂直切变为主的风切变区域，对飞行的危害程度较之强对流天气略弱(孙一昕和方娟，2012)。Williams(2005)提到辐射逆温型的低空急流天气可产生逆温风切变，由于其强度通常较小，容易被忽视，但若处置不当也会发生飞行危险;机场周围地物、建筑等也会引起风切变现象，且强度较大，这使机场的选址尤为重要(秦娟等，2011;迟继峰等，2012)。因此及时准确地判断风切变的位置、类型和强度，可有效减轻和避免危害，确保飞行及起降安全。

Chan and Shun(2010)研究认为风切变的强度标准涉及气象条件、飞机性能及飞行员驾驶技术等多方面因素，因此以其对飞行的危害程度作为强度分类标准成为航空界共识。按照国际民航组织建议(章澄昌，2000)，水平风的垂直切变强度标准见表1，其中空气层垂直厚度为30 m，时间尺度为2 min。当切变值达到0.1(1/s)以上时就威胁飞机安全。水平风的水平切变标准一般依据美国机场低空风切变报警系统所采用的报警标准(Britt et al．，1993)，水平风水平切变值在2.6(m/s)/km以上时极有可能对飞行构成危害。

表1 水平风垂直切变强度的标准Table 1 The standard ofverticalshear strength for horizontal wind

多普勒天气雷达不仅具有多普勒径向速度信息，还有反映风场变化的速度谱宽信息。在多普勒雷达资料的各种算法产品中，合成切变产品可反映大气流场的脉动性，同时可以展示3～10 km尺度以内危险的阵风锋、低空切变线等低空风切变现象以及与中尺度气旋有关的风切变识别产品(张京英等，2009;梁建宇和孙建华，2012)。在风切变识别方面，国内外有许多研究进展。目前主要的识别算法有直接计算差值滤波法、最小二乘法等(魏耀和张兴敢，2010)。Harris et al．(1985)首次提出了直接计算差值滤波合成切变算法。该算法在笛卡尔坐标系中分别计算一维径向风切变和一维方位风切变，然后经过求平均及滤波等步骤，合成为二维风切变。合成切变产品用于监测阵风锋、微下击暴流以及冷锋过境等风切变，并经逐步改进与机载雷达相结合，建立了早期的机载低空风切变预警系统。METSTAR公司生产的CINRAD-SA/SB等多普勒雷达均采用基于美国NEXTRAD的合成切变算法，但国内版权受限。在参考WSD-88D直接计算差值滤波合成切变算法的基础上，王珊珊等(2008)进行了雷达回波自动识别锋线的方法研究，产品能够直观识别风切变。胡明宝等(2000)用机场多普勒雷达资料根据最小二乘法计算合成切变值，定量描述机场附近的低空风切变区域，确定飞机起降的潜在危险区。2008年北京奥运会期间，气象保障部门应用最小二乘法获得的合成切变产品为人工影响天气作业提供了重要参考。魏耀和张兴敢(2010)对直接计算差值滤波和最小二乘法作对比分析，后者在定位精度、识别能力及边缘识别能力等方面具有优势。

机载多普勒雷达的风切变识别算法需要具有快速、简捷和高效的特点，本文根据机载气象雷达的特点和应用需求，采用最小二乘法的识别原理，通过实例分析，将此算法与目前多普勒天气雷达业务应用系统PUP上的风切变识别产品进行对比分析，讨论它在机载雷达上应用的可行性，为改进机载雷达风切变预警系统提供参考。

1 最小二乘法的识别算法原理

1.1 一维径向风切变

为了减少雷达测量噪声和奇异点的影响，选取合适的拟合窗口，其中包含n个距离库(同一径向上)，即(v1，r1)…(vi，ri)…(vn，rn)(vi为窗口中第 i个点的径向速度，ri为第i个点与雷达之间的距离)，然后利用最小二乘法计算雷达径向速度沿雷达径向的变化值CRs。设v*是拟合窗口的速度估计值，则线性回归方程为v*=α+βr。其中α、β是待定系数，α为初始值，β为斜率。

估计值与实测值之间的误差为

根据最小二乘法原理，使上述误差平方和最小，即令E最小，

则

解(3)式得

如图1所示，A、C两点之间的径向切变为(vC－vA)与资料点总长度ΔR之比，而斜率β是(vC－vA)与A、C两点之间距离Δr的比值。则计算一维径向切变值:

其中:LR是资料点采样体积在径向上的单位长度。

图1 径向切变示意图Fig．1 The illustration of radial shear

1.2 一维方位风切变

对于一维方位风切变(CAs)，在同一探测距离圆周上选取拟合窗口，计算过程中采取的拟合窗口内包含 n 个距离库，即(v1，θ1)…(vi，θi)…(vn，θn)(vi为窗口内第i个点的径向速度，θi为相应点的方位角)。由最小二乘法得拟合曲线的斜率β:

距离库内的数据是以该点为中心的采样体积内的平均值，采样体积与波束宽度和距离库长相关。如图2所示，同一库距内A、B两点径向速度的一维方位风切变值CAs(雷达的角度分辨率为1.5°)为:

图2 方位切变示意图Fig．2 The illustration of azimuth shear

由最小二乘法计算得 β=(vB－vA)/(rΔθ)，一维方位风切变值CAs为:

1.3 二维合成风切变

雷达径向速度沿径向、方位的综合变化，即一维径向风切变与一维方位风切变的合成，得到二维合成切变:

2 雷达资料处理

2.1资料质量控制

为保证识别效果，计算前需对雷达资料进行质量控制，包括雷达缺失数据填补及去除孤立噪声点等。缺失数据主要是雷达扫描时数据重置和原始数据丢失所致(Lakshmanan et al．，2003;袁杲等，2007)。为填补缺失径向数据，先判断方位角索引，若与上次相同则另外保存。读完后再判断是否有突变，若没有则不需进一步处理，否则就再对原始回波进行修正。去除孤立噪声时，以当前回波点为中心，取5×5的模板，计算无效反射率因子值所占百分比P，若P大于75%，则判断该点为孤立回波点，需去除(黄铃光，2011)。

2.2 拟合窗口选取

选取拟合窗口(表2)时，要兼顾随机载雷达探测距离远近而变化的径向切变和方位切变，以最大程度减少线性拟合的误差，因此为使合成切变的拟合精度不随距离而变化，规定组合拟合窗口面积保持不变。结合实验分析，资料点的选取应满足:

式中:Nr为径向切变拟合窗口资料点;Nt为方位切变拟合窗口资料点;ε为常数。

表2 径向切变、方位切变拟合窗口资料点数的选取Table 2 The data points of fitting window for radial shear and azimuth shear

鉴于机载雷达扫描方式、扫描半径与地基雷达的差别，在处理机载雷达资料时，拟合窗口的选取需要根据实际参数调整。

3 实例分析

3.1 石家庄机场风切变

1)个例1:2010年6月16—17日

2010年6月16日傍晚至17日下半夜，石家庄机场本场及周边地区出现大范围雷雨、短时暴雨、大风和冰雹等强对流天气，石家庄、沧州以北地区28个县(市)出现大风，万全地区瞬时风力最大达24 m/s;河北省除邢台东部、邯郸东部外均出现雷雨天气，中北部地区有33个县市降水量达25～106 mm;有8个测站出现冰雹，宽城测站冰雹最大直径为14 mm。这也是当年出现的范围最大、强度最强的对流天气过程。

图3是2010年6月16日20时(北京时，下同)的850 hPa天气图。可见，石家庄地区处于北部低压中心的东南侧，位于槽前脊后，降水量在20时前后并不多，但风速较大，东北西南走向的切变较明显。由于石家庄没有探空站，只能参考周边山西太原(53772)和河北邢台(53798)的探空资料。图4是2010年6月16日20时太原和邢台的温度对数压力(T-lgp)图，可见两站的对流有效位能CAPE值并不大，且水汽不充足，但此时850、700、500 hPa的风速较大，太原依次为 10、8、13 m/s，邢台依次为11、6、8 m/s，且邢台的 925 hPa 风速达 19 m/s，由此推断石家庄地区风速较大，风切变明显。

图3 6月16日20时850 hPa天气图Fig．3 850 hPa synoptic chart at 20:00 BST 16 June

图5为石家庄2010年6月16日23:54多普勒雷达径向速度资料图及PUP合成切变识别图和最小二乘法合成切变识别图。可见，图5b圆圈区域识别得模糊，而图5c相应区域则识别得清晰连续。

2)个例2:2010年7月10日

图46月16日20时太原(a)和邢台(b)的T-lgp图Fig．4 T-lgp plots of(a)Taiyuan and(b)Xingtai at 20:00 BST 16 June

图5 6月16日23:54雷达径向速度及合成切变图(半径:115 km;仰角:1.5°) a．PPI的雷达径向速度;b．PUP合成切变图;c．最小二乘法合成切变图Fig．5 Radar radial velocity and composite shear at 23:54 BST 16 June(radius:115 km;elevation:1.5°) a．radial velocity on PPI;b．composite shear of PUP product;c．composite shear of least square method

2010年7月10日石家庄及周边地区出现大范围雷雨、短时暴雨、大风和冰雹等强对流天气，石家庄机场本场20时左右遭遇大风天气，持续近半小时，最强阵风达20 m/s，降雨强、风速大给当地经济和人员造成重大损失。

图6是2010年7月10日08时850 hPa天气图。可见，河北省处于低压中心东南侧，西北干冷空气与东南暖湿空气在石家庄附近相遇并出现降水。20时低槽略东移，石家庄处于槽前脊后，降水增加，强对流明显(图略)。图7是20时太原和邢台的T-lgp图。可见，太原的对流有效位能 CAPE为658.3 J/kg，对流抑制能CIN为149.1 J/kg，邢台的对流有效位能 CAPE为256.9 J/kg，对流抑制能CIN 为39.4 J/kg。850、700、500 hPa的风速依次为:太原 4、7、16 m/s，邢台 1、10、15 m/s。分析表明石家庄将会继续存在风切变、降水和强对流天气。

图8为石家庄2010年7月10日21:24多普勒雷达径向速度资料图及PUP合成切变识别图和最小二乘法合成切变图。可见，图8b圆圈所示的弱切变区识别得较零散，而图8c则识别得较连续。

图6 7月10日08时850 hPa天气图Fig．6 850 hPa synoptic chart at 08:00 BST 10 July

图7 2010年7月10日20时太原(a)和邢台(b)的T-lgp图Fig．7 T-lgp plots of(a)Taiyuan and(b)Xingtai at 20:00 BST 10 July 2010

3)个例3:2009年6月5日

2009年6月5日午后安徽省大部出现雷雨大风及冰雹等强对流天气，至晚间上半夜全省40余县市出现8级以上大风，其中淮南20时53分最大瞬时风速达36 m/s，阜阳地区监测到强烈的阵风锋，怀远、淮南分别出现8、10 mm直径的冰雹。强对流引起大风天气在当地造成生命财产的严重损失。

图8 7月10日21:24雷达径向速度及合成切变图(半径:115 km，仰角:1.5°) a．PPI的雷达径向速度;b．PUP合成切变图;c．最小二乘法合成切变图Fig．8 Radar radial velocity and composite shear at 21:24 BST 10 July(radius:115 km;elevation:1.5°) a．radial velocity on PPI;b．composite shear of PUP product;c．composite shear of least square method

图9是2009年6月5日08时和20时的天气图。可见，西北干冷气流和西南暖湿气流在安徽交汇，该地区出现强对流并伴有雷暴天气。阜阳处于安徽北部，有西北—东南方向的强切变区，两时次天气图显示切变区稳定并维持，雷达速度图上切变带清晰可见。图10是20时阜阳T-lgp图。可见，图中红色区域所占比例较大，对流有效位能CAPE达1 782.7 J/kg;对流抑制能CIN为171.7 J/kg;温度露点仅为 3 ℃;925、850、700、500 hPa的风速依次为13、13、12、21 m/s，从地面到高空风速都很大，使得阜阳地区维持强对流天气和风切变。

图11为阜阳2009年6月5日21:44多普勒雷达径向速度资料图及PUP合成切变图和最小二乘法合成切变图。由图11b(圆圈处)和图11c均可见雷达西南侧存在阵风锋，而图11c识别的阵风锋边界比图11b更加清晰，说明最小二乘法识别的梯度强，边缘信息保留完整。

图10 6月5日20时阜阳T-lgp图Fig．10 T-lgp plot in Fuyang at 20:00 BST 5 June

图96月5日08时(a)和20时(b)850 hPa天气图Fig．9 850 hPa synoptic charts at(a)08:00 and(b)20:00 BST 5 June

3.2 识别结果分析

为进一步定量分析两种方法的识别效果，选取辐合辐散性气流(合成切变值点大于8.0×10－4s－1)在两种方法中对应的典型位置，进行定位误差分析。选取2010年6月16日23:54和2010年7月10日21:24的PUP合成切变图与最小二乘法合成切变图，分别进行对比统计，结果见表3和表4。选取2009年6月5日21:44的PUP合成切变图与最小二乘法合成切变图，对比统计阵风锋典型位置，结果见表5。

通过以上比较，得出如下结论:

图11 6月5日21:44雷达径向速度及合成切变图(半径:115 km，仰角:1.5°) a．PPI的雷达径向速度;b．PUP合成切变图;c．最小二乘法合成切变图Fig．11 Radar radial velocity and composite shear at 21:44 BST 5 June(radius:115 km;elevation:1.5°) a．radial velocity on PPI;b．composite shear of PUP product;c．composite shear of least square method

1)识别个数。对比分析PUP合成切变图和最小二乘法的合成切变图可知，在识别合成切变点个数时，最小二乘法识别得较多，且边缘切变的识别较突出，合成切变数值在8.0×10－4s－1以上的显示数量较多，参考价值高，在方位和距离上与速度图的辐散辐合位置较吻合;PUP产品显示的合成切变值个数较少且都在大值区，忽略了边缘切变的识别，参考价值不高。对比表明最小二乘法的识别能力优于PUP产品。

表3 2010年6月16日23:54两种方法识别结果的辐合辐散定位及误差分析(P为速度图中辐合辐散点的位置，P1为PUP产品中对应的位置，P2为最小二乘法对应的位置，D1为PUP产品定位误差，D2为最小二乘法对应误差)Table 3 Error analysis and location for convergence and divergence of two kinds results at 23:54 BST 16 June 2010(P is the location of convergence and divergence in radial velocity;P1 is the location of PUP product;P2 is the location of least square method;D1 is the location error of the PUP product;D2 is the error of least square method)

表4 2010年7月10日21:24两种方法识别结果的辐合辐散定位及误差分析(P为速度图中辐合辐散点的位置，P1为PUP产品中对应的位置，P2为最小二乘法对应的位置，D1为PUP产品定位误差，D2为最小二乘法对应误差)Table 4 Error analysis and location for convergence and divergence of two kinds results at 21:24 BST 10 July 2010(P is the location of convergence and divergence in radial velocity;P1 is the location of PUP product;P2 is the location of least square method;D1 is the location error of the PUP product;D2 is the error of least square method)

表5 2009年6月5日21:44两种方法识别结果的辐合辐散定位及误差分析(P为速度图中辐合辐散点的位置;P1为PUP产品中对应的位置;P2为最小二乘法对应的位置;D1为PUP产品定位误差;D2为最小二乘法对应误差)Table 5 Error analysis and location for convergence and divergence of two kinds results at 21:44 BST 5 June 2009(P is the location of convergence and divergence in radial velocity;P1 is the location of PUP product;P2 is the location of least square method;D1 is the location error of the PUP product;D2 is the error of least square method)

2)定位精度。表3—5表明，最小二乘法的定位与辐合辐散等合成切变大值点的实际位置差别较小，PUP产品定位误差较大。表3中D1值最大为0.8°、0.9 km，D2 值相应为 0.3°、0.1 km。表4 中，D1 值最大为 1.7°、0.5 km，D2 值相应为 0.1°、0.1 km;D1值最大为2.0°、0.8 km，D2值相应为0.4°、0.1 km。在定位精度上，最小二乘法较准确。

3)边缘信息。由对比分析看出，PUP合成切变产品的边缘信息不完整，边缘点的切变识别较少或直接过滤掉。最小二乘法的合成切变拟合效果较好，边缘信息提取完整，切变识别个数较多。分析2009年6月5日21:44的PUP合成切变图和最小二乘法切变图，可直观地看出PUP合成切变图中阵风锋较零散，出现较明显的锯齿状，而最小二乘法合成切变图上识别的阵风锋边缘清晰。

4)连续性。分析2009年6月5日阜阳资料，最小二乘法拟合特性保证了信息的完整，阵风锋的位置、平滑度、连续性等方面明显好于PUP合成切变产品，在分析对比连续时次的径向速度时，能够找到速度资料变化的依据，这为预报未来时次的危险切变区域提供了参考，有助于改进机载雷达风切变预警。

4 结论

为探索机载气象雷达的风切变识别算法，本文用地基多普勒雷达速度资料研究了风切变识别的最小二乘法，通过分析石家庄两次大风过程及阜阳一次阵风锋的实例，对比了PUP合成切变产品与最小二乘法合成切变的识别效果，结论如下:

1)在对多普勒雷达速度资料进行识别之前，结合相应时次的天气图和T-lgp图对天气形势做出预测非常重要。天气形势是产生局地强对流的背景条件，可大致确定对流性天气现象的位置，再结合雷达速度资料进行识别分析，有助于对风切变危险区域做出预警。

2)最小二乘法拟合效果较好，边缘信息较完整，径向速度资料的一维方位切变和一维径向切变识别与径向速度图中的切变、辐合、辐散信息相吻合。其二维合成切变在辐合或辐散等现象出现时，都会产生大值区域，这是对一维切变风场的必要补充。

3)三次个例分析表明，最小二乘法在识别效果、定位精度、边缘识别等方面均优于PUP合成切变产品。

4)从2009年6月5日阜阳多普勒雷达速度资料合成切变识别图上可以明显看出，在阵风锋的识别方面，最小二乘法具有良好的拟合特性，识别出的阵风锋在边界清晰度、平滑性及连续性等方面均具优势。

5)针对机载雷达资料的宽波束、风速等级少及向下扫描时远处距库高度低及易平滑的特点，探索机载雷达风切变算法时拟合窗口的选取尤为关键，需综合考虑机载雷达参数特点，以达到最佳识别和预警效果。

迟继峰，钟天宇，刘庆超，等．2012．复杂地形多测风塔综合地貌及风切变拟合修正的风资源评估方法研究［J］．华电技术，34(11):75-78．

胡明宝，谈曙青，汤达章，等．2000．单部多普勒天气雷达探测低空风切变方法［J］．南京气象学院学报，23(1):113-119．

黄铃光．2011．基于雷达数据的强对流天气的识别算法及实现［D］．南京:南京信息工程大学．

黄仪方．2002．航空气象［M］．成都:西南交通大学出版社:131-138．

梁建宇，孙建华．2012.2009年6月一次飑线过程灾害性大风的形成机制［J］．大气科学，36(2):316-336．

秦娟，吴仁彪，苏志刚，等．2011．机载前视气象雷达地杂波仿真方法［J］．现代雷达，33(8):72-75．

孙一昕，方娟．2012.2010年5月6日重庆强对流过程的天气学分析［J］．气象科学，32(6):609-621．

王丛梅，景华，王福侠，等．2011．一次强烈雹暴的多普勒天气雷达资料分析［J］．气象科学，31(5):659-665．

王珊珊，黄兴友，苏磊，等．2008．利用雷达回波自动识别锋线的方法［J］．南京气象学院学报，31(4):563-573．

魏耀，张兴敢．2010．多普勒天气雷达合成切变算法及改进方法的研究［J］．电子与信息学报，32(1):43-48．

武虎子，耿建中，罗振谊，等．2013．大型飞机在风切变环境中飞行安全尺度分析［J］．飞行力学，31(4):305-308．

袁杲，张志强，谢明元．2007．一种CINRAD雷达基数据读取丢失的改进方法［J］．成都信息工程学院学报，22(2):178-181．

张京英，孙成武，王庆华，等．2009．一次飑线大风的多种资料分析和临近预报［J］．气象科学，29(1):126-132．

张培昌，杜秉玉，戴铁丕．2001．雷达气象学［M］．北京:气象出版社:122-128．

章澄昌．2000．飞行气象学［M］．北京:气象出版社:179-187．

Britt C L，Harrah S D，Crittenden L H．1993．Microburst hazard detection performance of the NASA experimental windshear radar system［C］//AIAA93-3943:AIAA Aircraft Design，Systems and Operations Meeting．Monterey，CA:AIAA:1-13．

Chan P W，Shun C M．2010．Latest developments of wind-shear alerting services at the Hong Kong International Airport［C］//14th Conference on Aviation，Range，and Aerospace Meteorology．Atlanta，Georgia:American Meteorological Society:17-21．

Harris F I，Glover K M，Smythe G R．1985．Gust front detection and prediction［C］//Preprints，14th Conference on Severe Local Storms．Boston:Bulletin of the American Meteorological Society:342-345．

Lakshmanan V，Hondl K，Tumpf G，et al．2003．Quality control of WSR-88D data［C］//Preprints，31st International Conference on Radar Meterology．Seattle，WA:American Meteorological Society:522-525．

Williams P．2005．Aircraft trajectory planning for terrain following incorporating actuator constraints［J］．Journal Aircraft，42(5):1358-1362．