一种大气风场急流线计算方法

甘建红， 漆 慧,胡文东， 舒红平， 罗 飞， 何童丽,尹麒名， 黎仁国

(1.成都信息工程大学， 成都 610200； 2.西华师范大学， 南充 637002)

1 引 言

大气风场急流现象最早在20世纪30年代被人们发现，急流与对强对流天气形成发展、风能转化、空气污染物的输送和扩散、森林火灾的蔓延等有着密切联系[1-2].目前实际气象业务中通过人工手动方式在软件上进行识别和绘制急流，一定程度上限制了过程自动化、预报智能化的快速发展.随着智能天气预报和格点预报兴起，计算机自动识别并标识气象风场急流的研究逐渐受到人们重视.

Molnos等[3]以风速、风向和纬度导向函数构造代价函数，通过Dijkstra方法寻找最短路径方式计算副热带急流(Subtropical Jet Stream)和极锋急流(Polar Jet Steam)，由于该急流线是两条的闭合线，没有考虑分叉和急流合流的现象，且算法需要已知起点，可见该算法不适合于非全球范围风场数据中计算具有多个急流的情况.

金宏忆[4]以VAD算法为基础，用雷达探测资料以及多普勒天气雷达二次产品，探讨了机器识别低空急流的检测方法.章晗等人[5]也利用雷达资料研究急流识别，该算法只针对西南风低空急流自动识别做出了研究，未对其它类型的急流提取回波特征.

路志英等[6]和王萍等[7]的方法类似，对站点风数据计算急流线并用MICAPS进行显示，该方法利用地面站的观测资料，经过预处理、聚类及拟合实现了850 hPa 等压面上中国东南部沿海地区低空急流的绘制.该方法根据该地区以东南风和西南风为主的特征进行识别，对其它地方适用性需要进一步优化.

随着智能网格预报的兴起，不但预报数据采用格点形式，实况数据也逐渐采用格点数据形式进行应用[8-9].本文针对格点类型的风场数据进行急流轴计算，首先通过方向相干性计算急流区域，然后采用骨架提取算法获取接近几何中心的格点，基于获取格点位置采用重心公式计算急流轴的点，最后连接急流轴的点得到急流轴.

2 急流计算方法

算法主要经过三个阶段实现：计算急流区域，计算中心轴(骨架)点及其微调，连接点得到急流线，具体算法流程如图1所示.

3 计算急流区域

大气风场急流具有风速较大、局部区域内风向具有一致性特点.对于站点数据可采用图的方法计算点间的相似性，而本文输入格点数据为棋盘格分布的数据，因此采用Hessian的方法计算相邻点的相干系数.本文分两步计算激流区域：(1) 计算风向一致性好的区域；(2) 从得到区域中去掉风速小的地方，保留下来的区域即为急流区域.

为了便于使用算子进行计算，在格点数据的上边增加一行数据，其值等于最边上一行数据，同理在格点的下边增加一行数据，在左边、右边也各增加一列数据，将增加数据之前的格点称为原格点.

图1 急流自动识别流程图Fig.1 Jet automatic identification

(1) 计算风向一致性较好的区域.

首先，对全部原格点数据计算方向相干性系数.在图像处理等领域内相干性系数计算主要面向标量场，而风场是一个矢量场，为了能让其刻画方向保持一致性特性，提出算子范围内水平、竖直方向的最大风向夹角作为X、Y方向差分并构建Hessian矩阵[10]的方法计算相干性系数，得到Hessian矩阵在此表达式为

(1)

其中，Fx=max{θ|邻域内水平方向两两格点风向夹角};Fy=max{θ|邻域内竖直方向两两格点风向夹角}.

具体离散点中计算方法如下：假设当前计算风向相干系数格点位置为(x0,y0)，本文实验中邻域尺寸为3×3，如图2所示.计算(x0,y0)点风向与(x0+1,y0)点风向构成的夹角，记为Fx1；计算(x0,y0)点风向与(x0-1,y0)点风向构成的夹角，记为Fx2；计算(x0+1,y0)点风向与(x0-1,y0)点风向构成的夹角，记为Fx3.于是,Fx1、Fx2、Fx3的最大值即为Fx.同理，得到y方向的差分值Fy.

然后，计算该Hessian矩阵的特征值λ1和λ2，并按照式(2)计算风向相干性系数k.如果原格点的系数小于阈值T，则为急流区域候选点.

k=|λ1-λ2|

(2)

(2) 提取急流区域.

按照公式(3)计算每个候选格点风矢量的模值，即风大小，如果模值大于阈值S，则为急流区域内的格点，否则不是急流区域的格点.至此可得到一个或多个急流区域.

(3)

其中,vx、vy分别为风矢量在水平、竖直方向的分量.

(a) (b) (c)

4 计算急流区域几何中心轴

《中尺度天气图技术规范(修改稿)》，低空急流的技术要求是：当有两个以上连续测站风速超过12 m/s时， 沿12 m/s以上大风区的几何中心分析低空急流.基于该规范的思想，算法在计算出急流区域后利用骨架提取方法计算区域的骨架点，骨架点即也是格点.

4.1 计算急流轴附近的格点

骨架提取算法[11-12]是提取几何中心轴常用的方法，而查表法是骨架提取的一个重要方法.查表法将一系列约束条件用周围像素权重乘积之和通过查表方式进行表达，比如内部点不能删除、孤立点不能删除、直线端点不能删除等约束条件.在此，每个格点的值只能为1或0，如果为急流区域的点，则为1，否则为0.权重分布如图3所示，对于任意一格点，取其周围紧邻8个格点为邻域，并按照一定顺序给与权重，权重分别为1、2、4、8、16、32、64、128，中间方格表示当前格点.

图3 邻域权重分布图

计算方法为将当前格点周围8个格点的值分别乘以对应位置的权重，并对乘积之和除以255得到最终结果，在表1中查询该最终结果对应的值，对应值如果是1，则表示这个点应该删掉，反之若值为0，表示这个点应该保留.

4.2 微调急流轴点

通过上述算法获取的急流轴位于格点之上，但在实际情况中，风场的急流轴刚好经过格点的概率非常小.采用重心公式对骨架点进行微调，使得骨架点接近风速最大位置，与实际情况一致.具体方法为，设第i个骨架点为当前待微调的骨架点，将当前骨架点及其左、右、上、下、左上、右上、左下、右下共9个格点的风矢量的模带入式(4)和(5)，9个格点中如果某个格点不在急流区域将其风大小F设置为0进行计算，得到微调后的骨架点位置(Xi，Yi).

(4)

(5)

其中,xij，yij分别是当前骨架点的横坐标与纵坐标;Fij是格点的风矢量的模.图4(a)中圆点标注的是计算得出的急流区域，五角星标注的是微调后的骨架点.

最后,将计算的急流轴点(微调后的骨架点)进行连接得到急流线，如图4所示，通常需要对其进行平滑处理，平滑处理不是本文研究重点且方法较为成熟，不在赘述.

(a)

(b)

5 算法测试结果及分析

算法实验采用MeteoInfo开源组件读取数据，在VisualStudio2015开发环境中使用C#语言实现，输入数据为MICAPS第11类数据，程序结果为急流轴的经纬度坐标.

5.1 统计分析

将算法应用于2010年1月1日至7月1日及2017年7月1日至12月31日08～20时500 hPa风场数据中的352组数据进行测试.将测试结果交若干位专家进行核实，如果专家认为与手动绘制较吻合则定义为准确击中；如果与人工绘制偏差较大但是流线方向与风向趋势一致则定义为偏差击中；若本算法识别存在低空急流但人工分析不存在低空急流则认为空报；反之，若人工分析存在低空急流但本算法并未识别出，则认为漏报.实验结果如表2所示.

图5是2010年1月1日至4日08时高度为500 hPa的UV数据实验结果，红色线条为计算的急流线，蓝色为表示风向的风廓线.算法获得的急流走向和风向走势具有良好的一致性.

表2 算法应用于2012年和2017年部分数据结果

(a)

(b)

(c)

(d)

5.2 个例分析

将本算法应用于2017年7月14日8时150 hPa的UV数据,结果如图6所示.图b中圆形点为急流区域内的格点，圆形格点组成的区域即为急流区域.五角星为微调后得到的骨架点，粗线即为急流轴.通过自动检测能快速检测出急流轴，对于小短的急流轴也能找到，而这些通常是人工识别容易忽略的急流轴.急流的合并与分流是描述风场变化的重要特性，算法计算的急流轴也保持了合并与分流的情况.

文献[7]对2009年3月21日8时、2009年4月19日8时地面站数据进行了急流线分析，本文针对2009年3月21日8时、2009年4月19日8时的实况格点数据进行实验，本文算法的急流线在北海附近更靠近最大风速、更靠近风力大的站点，与急流线的定义更一致.由于文献[7]基于站点数据计算急流，因此本文算法急流经过海洋较多，文献[7]的急流主要在陆地，由于数据源的不同，结果存在一定差异.

图6 高空风场计算结果Fig.6 Result of high level UV

一定程度上，风场本身也是流体，根据流线的性质，二维平面内同一时刻的不同流线不能相交,因此，本实验结果没有出现急流相交的情况，更具有可解释性，符合流体力学性质.

6 结 论

本文根据大气风场格点矢量数据，利用相干性系数获取特征值，将格点矢量数据转为标量.再利用骨架提取的方法可以自动识别急流区域并绘制急流线.相比人工处理效率更高、准确度更高，为天气预报自动化、智能化提供方法参考.为了得到更好的结果，可对原始格点数据进行插值预处理，比如基于径向基的插值方法.

根据《中尺度天气图技术规范(修改稿)》，低空急流的技术要求是：当有两个以上连续测站风速超过12 m/s时， 沿12 m/s以上大风区的几何中心分析低空急流.可见算法首先计算急流区域然后通过骨架提取几何中心轴与该规范思路一致.