CN106919792A - 基于高精度数值风场资料的涡旋中心自动识别方法 - Google Patents

Abstract

本发明公开了一种基于高精度数值风场资料的涡旋中心自动识别方法，包括：1)依据850pha数值风场的矢量特征，利用衡量风向杂乱度的特征提取出包含轴对称型风场的感兴趣区域；2)利用圆形数据理论对理想的轴对称型风场进行特征描述和分类；3)相对于复杂的轴对称风场依据步骤2)的结论和构建变形系数，初步检索出涡旋候选中心点集；4)改进Ward聚类算法得到属于同一涡旋系统的中心点集；5)利用步骤3)的候选中心点集的涡度对涡旋系统分类，并精确定位全局的涡旋中心。本发明可自动识别和定位多种类型的轴对称风场中心，能识别多尺度的涡旋系统及复杂风场中带有一定形变的涡旋系统；从而实现涡旋系统的准确、完整自动识别和涡旋中心定位。

Description

基于高精度数值风场资料的涡旋中心自动识别方法

技术领域

本发明涉及气象学领域，特别涉及一种基于高精度数值风场资料的涡旋中心自动识别方法。

背景技术

在中尺度天气系统分析领域，基于地面和高空的风场数据是预报员在分析和跟踪天气系统的主要资料依据。其中，辐合(或辐散)式涡旋是构成强对流或灾害性天气的重要风场系统。例如在北半球，中尺度的逆时针涡旋是灾害性风暴,龙卷风的重要成因；大尺度的热带气旋，是形成我国夏季台风的主要因素；反气旋也是形成雷暴天气的重要载体。因此，涡旋自动识别算法对于上述灾害性天气的预警和预测具有至关重要的作用(Corey K.和Potvin 2013)[1]；并且，实现计算机对涡旋的自动识别和量化描述，会有助于智能处理目前地面和高空庞大的数值资料，有助于利用数据挖掘技术在大数据中探索新规律。

目前，国内外气象界研制开发的中小尺度的涡旋检测算法大部分是基于雷达数据的。(Potvin等2009[2]，Potvin等2011[3])。对于中大尺度的涡旋识别，大都基于网格气压数据的局部高(低)气压中心这一判别条件(John和Rodrigo 2012[4]，林志强，周振波2013[5])。基于网格数据的风场涡旋识别算法在气象领域相对较少：Mark,R.和Sinclair.(1994)提出基于风场的位势相对涡度代替气压低中心[6]。但是，高涡度值的风场包含了涡旋区，风切变区以及风向变化杂乱的区域(Corey K.和Potvin 2013)[1]，因此单纯地计算高涡度区，对于检测气旋会产生大量的空报。Mark和Sinclair(1997)[7]提出了基于上述方法的扩展算法：通过一种空间平滑技术来控制需检测涡旋的尺度，保证了在不同分辨率的网格数据下对固定尺度的气旋检测的一致性。同时，结合气旋的中心低气压值和高涡度值的特性作为衡量气旋强度的条件。但其仍然没有解决空报的问题。Naylor和Gilmore(2012)[8]通过在网格风场数据上添加规则(包括：垂直涡度，水平风速的大小和中心点向四周方向的水平气压梯度)识别龙卷风。但该方法疏忽了实际的复杂流场上涡旋不具有完整的轴对称结构的现实，也就是相对于理想的圆形涡旋结构，实际的涡旋往往带有不同程度的形变。针对此问题，Corey K.和Potvin(2013)[1]提出针对形变涡旋结构的对流性涡旋的检测模型：将基于多普勒雷达数据的风暴检测和对流性涡旋检测方法应用在笛卡尔网格风场数据中。该模型由兰金涡旋模式和恒定流场、线性切变流场和线性辐散流场线性组合而成。然而，受气压、温度、地势等多环境因素的影响，实际的风场在中大尺度下已不具有线性的性质，当流场较复杂时，模型的误差较大，只有缩小检索域尺度才能降低误差。

[参考文献]

[1]Corey K.Potvin.A Variational Method for Detecting andCharacterizing Convective Vortices in Cartesian Wind Fields.Monthly WeatherReview.2013,141(9):3102-3114。

[2]Potvin,C.K.,A.Shapiro,T.-Y.Yu,J.Gao,and M.Xue,Using a low-ordermodel to detect and characterize tornadoes in multiple-Doppler radardata.Mon.Wea.Rev.,2009,137:1230–1249。

[3]Potvin,C.K.,A.Shapiro,M.I.Biggerstaff,and J.M.Wurman,The VDACtechnique:A variational method for detecting and character-izing convectivevortices inmultiple-Doppler radar data.Mon.Wea.Rev.,2011,139:2593–2613。

[4]John Hanley,RodrigoCaballero.Objective identification and trackingof multicentre cyclones in the ERA-Interim reanalysis dataset.QuarterlyJournal of the Royal Meteorological Society.2012:612-625。

[5]林志强，周振波假拉.高原低涡客观识别方法及初步应用.高原气象[J],2013,32(6)。

[6]Mark,R.Sinclair.An Objective Cyclone Climatology for the SouthernHemisphere.Monthly Weather Review[J].1994:2239-2256。

[7]MARK,R.S INCLAIR.Objective Identification of Cyclones and TheirCirculation Intensity,and Climatology.WEATHER AND FORECASTING.1997:595-612。

[8]Naylor,J.,and M.S.Gilmore,Environmental factors influential to theduration and intensity of tornadoes in simulated supercells.Geophys.Res.Lett.,2012,39。

发明内容

本发明公开了一种基于高精度数值风场资料的涡旋中心自动识别方法，可以解决的技术问题包括：自动识别和定位多种类型的轴对称风场中心(包括：顺(逆)时针旋转涡旋，辐合(散))；能识别多尺度的涡旋系统；能够识别复杂风场中带有一定形变的涡旋系统；达到准确、完整自动识别涡旋系统和定位涡旋中心的目的。

为了解决上述技术问题，本发明提出的一种基于高精度数值风场资料的涡旋中心自动识别方法，包括以下步骤：

步骤一、依据850pha数值风场的矢量特征，利用衡量风向杂乱度的特征提取出包含轴对称型风场的感兴趣区域；步骤如下：

1-1)设输入的风场矩阵为W_M,N，其中，M,N分别为矩阵的行数和列数；设风场矩阵W_M,N中的两个矢量的方向夹角为Λ，构建一个统计特征即：在L×L窗口内，分别计算风场矩阵W_M,N中的网格点(i,j)和该网格点的相邻点(i，j+1)与(i-1,j)的方向夹角Λ，并累积求和；设风矢量的角度分量为θ，则两个矢量的方向夹角为Λ(θ₁,θ₂)，统计特征的计算公式如下：

式(1)中，i，j，m和n分别是风场矩阵W_M,N中网格的东西方向和南北方向的索引；m＝r+2，…，M-r；n＝r+1，…，N-(r+1)，L为窗口的尺度，且为奇数；

1-2)对上述统计特征进行归一化得到归一化的统计特征对输入的给定分辨率的高精度网格风场数据，计算出统计特征的最大值，即对进行归一化得到计算公式为最后，遍历计算W_M,N，得到归一化的特征矩阵即为：

1-3)利用Bayesian决策对参数L进行优化，和对公式(1)得出的进行阈值选取，以获得输入风场矩阵的最优感兴趣区；从给定分辨率的风场数据中搜集涡旋样本V，等分为三组测试样本，包括：V_A，V_B和V_C，其中，风场数据包含不同季节、时间和海拔高度；然后，对杂乱风样本M的三个子样本：M_A，M_B和M_C和一致风样本S的三个子样本：S_A，S_B和S_C分别计算当L＝5,7,9,11和13的五类训练子特征：和当L取上述给定某一尺度值时，分别计算对于类别V_A，M_A和S_A的条件概率密度：和根据最小错误率的原则，获得的最优阈值和相应的最小错误率同样的过程应用于[V_B,M_B,S_B]和[V_C,M_C,S_C]，获得相应的最小错误率和比较不同的L取值下的最小错误率的平均值，求得对于给定分辨率的风场数据的最优L取值；基于最优L值，确定的最优阈值，该阈值包含一个与杂乱区域分界的最大值和与一致风向区域的最小值。

步骤二、利用圆形数据理论，对理想的轴对称型风场，包括：顺时针旋转场、逆时针旋转场、辐合场、辐散场和鞍形场进行特征描述和分类；步骤如下：

2-1)对于两个圆形数据变量α和β，它们之间的一种有向的方向差Q(α,β)，定义为：

定义α到β的一种走向为：

其中，α和β的取值范围各为[-180°,180°]；α到β的走向α→β为：在极坐标系下从α到β的走向；

2-2)结论1：设θ_k是以理想轴对称风场中心为中心的单位圆弧上k位置处风的角度，当0°<Q(θ_k,θ_k+ε)<180°,ε>0时，则θ_k→θ_k+ε:counterclockwise；

对于所有的0°≤k<360°，结论1都成立时，判定θ_k为第一类轴对称场，即涡旋场或辐合场或辐散场；

结论2：设θ_k是以理想轴对称风场中心为中心的单位圆弧上k位置处风的角度，当-180°<Q(θ_k,θ_k+ε)<0°,ε>0时，则θ_k→θ_k+ε:clockwise；

对于所有的0°≤k<360°，结论2都成立时，判定θ_k为第二类轴对称场，即鞍形场；

2-3)设是轴对称型风场，对于所有的0°≤k<360°，当时，和θ_k同属一类轴对称场；如果θ_k为逆时针涡旋场，当δ＝90°时，为辐合场；当δ＝180°时，为顺时针涡旋场；当δ＝270°时，为辐散场；如果θ_k为鞍形场，δ取0°～360°任意值，为鞍形场；

步骤三、相对于理想轴对称风场，实际的涡旋场是复杂的轴对称风场，依据步骤二的结论和构建变形系数，初步检索出涡旋候选中心点集；步骤如下：

3-1)对步骤一计算出的感兴趣区域的每个格点风向值进行角度量化：依次为：东风、东北风、北风、西北风、西风、西南风、南风和东南风八个方向；

3-2)以步骤一计算出的感兴趣区域为当前检索范围，依次对每个以当前格点为中心的d×d窗口下分析域的风向进行特征提取，判定当前格点是否为涡旋中心点；从中心点向X轴方向的射线开始，以45°角为步长，依次引出八条射线对分析域进行八等分；

3-3)涡旋场的特征提取；对于当前分析域的第i(i＝1,2,…,8.)个子区域，所有格点数为S_i、具有相同角度α的连通格点数为则该子区域的特征值定义为c_i，因此当前分析域的特征向量描述为C＝{c₁,c₂,…c₈}；C的提取规则为：

IF存在α使得

c_i＝α；

ELSE

停止并检索下一格点；

END

其中，ρ是阈值变量，默认值是0.4；

3-4)将C简化为设c_i和c_j是两个相邻子区域的特征值；即，其中，如果i<8,j＝i+1,如果i＝8,j＝1；简化规则为：

IFc_i＝c_j

去掉c_j

END

3-5)剔除不合理的设q是的长度，和是的相邻元素，其中，如果i<8,j＝i+1,如果i＝8,j＝1；则规则为：

IF q<4或存在c_i使得

停止并检索下一格点；

END

3-6)设为中方向角所对应的空间子象限的索引向量；对通过步骤3-5)的和映射为单位圆的内接多边形；规则是：以或中元素的值作为单位圆的内接多边形的顶点；或映射的19种圆内接多边形如表1所示；

表1

3-7)对步骤3-6)中的所有多边形构建形变指数λ，去除形变过度的候选中心，并获得确定下的涡旋的形变指数，如表1所示；λ的构建方法如下：

设圆形数据向量A＝{α₁,α₂,...,α_n}，其中，α_i≠α_j，当0°<Q(α_i,α_i+1)<180°且0°<Q(α_n,α₁)<180°时，构建形变指数λ；λ是圆形数据向量A的合成向量长度、圆形数据向量A均匀波动性和圆形数据向量A的空缺率的线性组合；

(1)合成向量长度：r_o

r_o的取值范围为[0，1]；r_o越接近1，α₁,α₂,...,α_n越聚焦于它们的均值；

(2)均匀波动性：δ_o

圆形数据向量A的有向方向差为Q(A),则Q(A)的标准差为：

同样地，表1中所有圆内接多边形的δ值中的最大值为1.73，得出标准化后的δ_o：

(3)空缺率：ρ_o

设圆形数据向量A的维度为n,空缺率ρ_o的计算如下：

最后，将上述三种特征线性组合，得到一个综合的形变指数λ：

λ＝ar_o+bδ_o+cρ_o (9)

式(9)中，0<a<1,0<b<1,0<c<1且a+b+c＝1.经实验得出a＝b＝0.25,c＝0.5，所有圆内接多边形的λ值如表1所示；

3-8)排除涡旋的伪中心点

从风场数据中挑选涡旋样本估计λ的值；首先，每个样本的和分别由d＝5,7,9的分析域计算得到；利用步骤1-3)中的Bayesian决策得到最优λ值，从而得到：

或

停止并检索下一格点；

END

步骤四、改进Ward聚类算法，对属于同一个涡旋系统的中心点集聚类；步骤如下：

4-1)对步骤三使用d＝5，7和9对感兴区域进行多尺度窗口的搜索候选中心点集P，其中P包含了多个涡旋系统的候选中心点集；利用Ward聚类算法将P中的各点聚类到各自的涡旋系统得到P_k，k＝1,2,…,K；K为类别数；

4-2)对于类别数K进行自适应选取；步骤4-1)中，每步类合并计算的最小离差平方和为D_i，计算D_i在最大离散曲率T_i处对应的类别数为最终认定的K；T_i的计算方法如下：

T_i＝1-|cosω_i| (10)

式(10)中，ω_i＝angle(p_i-1,p_i,p_i+1)是线段[p_i-1,p_i]和[p_i,p_i+1]的夹角，而K＝argmax_i{T_i}；

步骤五、利用步骤三的候选中心点集的涡度对涡旋系统分类，并精确定位全局的涡旋中心；

5-1)得到了每类涡旋系统的候选中心点集P_k后，采用数学形态学中的膨胀算法对P_k进行膨胀处理，以完全覆盖涡旋中心区域，方法如下：

式(11)中，B是一个3×3的结构元素；

5-2)将涡旋类别判定为顺时针方向或逆时针方向；引进二维连续风场的垂直涡度ζ：

将式(12)转换为离散风向场的垂直涡度

式(13)中，θ_i,j是单位矢量的角，i和j是水平方向和竖直方向的索引；令l₁(或l₂)记为Z_k元素的正或负涡度的个数；l₁或l₂定义为如下：

式(14)中，#{·}集合{·}的基数；如果l₁>l₂，当前第k个涡旋判定为逆时针涡旋，否则判定为顺时针涡旋；

5-3)根据涡度值确定涡旋中心参考点；在涡旋中心区域Z_k，由公式(13)计算每个格点的垂直涡度将垂直涡度绝对值最大的检测点作为涡旋中心的参考点；

5-4)根据零风速约束精确定位涡旋中心；在到涡旋中心参考点距离小于ε的区域内，(ε＝3)，将网格划分为两个三角网格，在三角网格内通过线性插值算法计算出零风速点坐标位置，连接两个零风速点即可得到零等值线段；由此得到风速东西分量U和南北分量V在一个网格内的零等值线段，同一网格上U、V的等值线段的交点，即为风速零点；若整个局部区域存在多个风速零点，则选择距离气旋中心参考点最近的点作为气旋中心；若整个局部区域不存在风速零点，则对风速分量U、V分别进行线性插值，求解风速矢量U、V的强度极小值，将气旋中心点校正到风速矢量U、V极小值点位置。

与现有技术相比，本发明的有益效果是：

气象中的涡旋本身就属于风场的一个基本特性，由于在风矢量场上识别涡旋存在一定难度，目前国内现有的技术基本都是基于低气压中心识别气旋与反气旋中心。本发明的涡旋识别算法完全基于高精度的风场网格数据，通过构建风向纹理特征，提取出包含识别目标的感兴趣区域；利用多尺度窗口，提取带有空间位置信息的子区域风向向量，在此基础上建立一系列规则，包括圆形数据的走向、风向向量的形变度，搜索出候选涡旋中心点集；再将其聚类到各自的涡旋系统；最后对涡旋的旋转模式分类，并定位出涡旋系统的中心点。本方法特别针对高精度的复杂风场网格资料。定位准确，并对带有形变的涡旋进行度量和筛选。最后，实验验证了本方法的有效性。

附图说明

图1(a)为使用L×L的窗口在风矢量矩阵W_M,N(实线矩形表示)上计算出特征矩阵Γ^L(虚线矩形表示)的过程示意。

图1(b)是一幅850pha的0.25°分辨率的风向场和灰度图的叠加图。颜色越深代表风向纹理越一致，颜色越浅代表风向纹理越杂乱。

图1(c)是(b)经过参数优化后的感兴趣区域(灰色区域)提取的效果图。其中，L＝9,

图2(a)到图2(e)是五类典型的中心对称型风场的简化图，其中：图2(a)逆时针旋转场的简化图，图2(b)是辐合场的简化图，图2(c)是顺时针旋转场的简化图,图2(d)是辐散场的简化图，图2(e)是鞍形场的简化图。

图3(a)是搜索到风场网格点(214，286)，并以该点为中心的一个9×9的分析域。

图3(b)将图3(a)的分析域八等分，在等分线上的网格点同属于相邻的两个子域。

图3(c)将图3(b)提取的特征向量映射到单位圆的内接多边形(实线标记)。

图3(d)将图3(b)提取的特征向量映射到单位圆的内界接多边形(实线标记)。

图4为是一幅850pha的0.25°分辨率的风向场通过d＝5,7和9窗口分析后，得到的候选中心点集。

图5是分别显示了逆时针涡旋和顺时针涡旋经过多尺度d的分析和规则验证，得到候选中心点集，以及通过方向涡度场对涡旋旋转方向的分类和全局中心点的定位。圆圈表示候选中心点，方框表示顺时针涡旋中心，三角表示逆时针涡旋中心。

图6是得到的涡旋中心定位点结果示意，该网格资料的高度为：850pha；经纬度范围为：60-150E，60N-10S；网格分辨率为：0.25°；图中D表示逆时针涡旋中心，G表示顺时针涡旋中心，d和g代表逆时针涡旋和顺时针涡旋中心参考点。

图7是本发明方法的流程图。

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施例仅对本发明进行解释说明，并不用以限制本发明。

本发明提供了一种基于高精度网格点风场资料的涡旋中心自动识别方法，本方法能自动检测出顺(逆)时针涡旋中心，对基于数值场的天气系统智能预报和统计分析，奠定了重要的研究基础。

实施例：给出一张包含我国主要地区的高度为850pha风场流线图，如图1(b)。其中，涡旋中心以人工标定，用来检验预处理的效果。

本发明基于高精度网格点风场资料的涡旋中心自动识别方法，如图7所示，包括以下步骤：

步骤一、依据850pha数值风场的矢量特征，利用衡量风向杂乱度的特征提取出包含轴对称型风场的感兴趣区域；步骤如下：

1-1)设输入的风场矩阵为W_M,N，其中，M,N分别为矩阵的行数和列数；设风场矩阵W_M,N中的两个矢量的方向夹角为Λ，构建一个统计特征即：在L×L窗口内，分别计算风场矩阵W_M,N中的网格点(i,j)和该网格点的相邻点(i，j+1)与(i-1,j)的方向夹角Λ，并累积求和；设风矢量的角度分量为θ，则两个矢量的方向夹角为Λ(θ₁,θ₂)，统计特征的计算公式如下：

式(1)中，i，j，m和n分别是风场矩阵W_M,N中网格的东西方向和南北方向的索引；m＝r+2，…，M-r；n＝r+1，…，N-(r+1)，L为窗口的尺度，且为奇数；遍历整个W_M,N，获得特征矩阵Γ^L的过程如图1(a)所示。

1-2)对上述统计特征进行归一化得到归一化的统计特征对输入的给定分辨率的高精度网格风场数据，计算出统计特征的最大值，即对进行归一化得到计算公式为最后，遍历计算W_M,N，得到归一化的特征矩阵即为：

1-3)利用Bayesian决策对参数L进行优化，和对公式(1)得出的进行阈值选取，以获得输入风场矩阵的最优感兴趣区；从0.25°分辨率的风场数据(包含不同季节、时间和海拔高度)中搜集涡旋样本V，等分为三组测试样本，包括：V_A，V_B和V_C；然后，分别计算当L＝5,7,9,11和13的五类训练子特征：和对杂乱风样本M的三个子样本：M_A，M_B和M_C和一致风样本S的三个子样本：S_A，S_B和S_C分别进行同样的计算和的操作。

当L取上述给定某一尺度值时，分别计算对于类别V_A，M_A和S_A的条件概率密度：和根据最小错误率的原则，获得的最优阈值和相应的最小错误率同样的过程应用于[V_B,M_B,S_B]和[V_C,M_C,S_C]，获得相应的最小错误率和比较不同的L取值下的最小错误率的平均值，求得对于0.25°分辨率的风场数据的最优L取值；经实验，求得对于0.25°分辨率的风场数据的最优L取值为9，最优阈值为图1(b)是0.25°分辨率的风矢量场和灰度图的叠加图，图1(c)是根据最优阈值求得的包含轴对称型风场的感兴趣趣区域。

步骤二、利用圆形数据理论，对理想的轴对称型风场(包括：顺时针旋转场、逆时针旋转场、辐合场、辐散场和鞍形场)进行特征描述和分类；步骤如下：

2-1)对于两个圆形数据变量α和β，它们之间的一种有向的方向差Q(α,β)，定义为：

定义α到β的一种走向为：

其中，α和β的取值范围各为[-180°,180°]；α到β的走向α→β为：在极坐标系下从α到β的走向；

2-2)结论1：设θ_k是以理想轴对称风场中心为中心的单位圆弧上k位置处风的角度，当0°<Q(θ_k,θ_k+ε)<180°,ε>0时，则θ_k→θ_k+ε:counterclockwise。

对于所有的0°≤k<360°，结论1都成立时，判定θ_k为第一类轴对称场，即涡旋场或辐合场或辐散场。

结论2：设θ_k是以理想轴对称风场中心为中心的单位圆弧上k位置处风的角度，当-180°<Q(θ_k,θ_k+ε)<0°,ε>0时，则θ_k→θ_k+ε:clockwise。

对于所有的0°≤k<360°，结论2都成立时，判定θ_k为第二类轴对称场，即鞍形场。

2-3)设是轴对称流场，对于所有的0°≤k<360°，当时，和θ_k同属一类轴对称场。如果θ_k为逆时针涡旋场，如图2(a)所示；当δ＝90°时，为辐合场，如图2(b)所示；当δ＝180°时，为顺时针涡旋场，如图2(c)所示；当δ＝270°时，为辐散场，如图2(d)所示；如果θ_k为鞍形场，如图2(e)所示，δ取0°～360°任意值，为鞍形场。

步骤三、相对于理想轴对称风场，实际的涡旋场是复杂的轴对称风场，依据步骤二的结论和构建变形系数，初步检索出涡旋候选中心点集；步骤如下：

3-1)对步骤一计算出的感兴趣区域的每个格点风向值进行角度量化：依次为：东风、东北风、北风、西北风、西风、西南风、南风和东南风八个方向。

3-2)以步骤一计算出的感兴趣区域为当前检索范围，依次对每个以当前格点为中心的d×d窗口下分析域的风向进行特征提取，判定当前格点是否为涡旋中心点，如图3(a)所示。从中心点向X轴方向的射线开始，以45°角为步长，依次引出八条射线对分析域进行八等分，如图3(b)所示。

3-3)涡旋场的特征提取：对于当前分析域的第i(i＝1,2,…,8.)个子区域，所有格点数为S_i，具有相同角度α的连通格点数为则该子区域的特征值定义为c_i，因此当前分析域的特征向量描述为C＝{c₁,c₂,…c₈}；C的提取规则为：

IF存在α使得

c_i＝α；

ELSE

停止并检索下一格点；

END

其中，ρ是阈值变量，默认值是0.4；对于图3(b)的实例，C＝{-90°,0°,0°,45°,135°,±180°,-135°,-90°}。

3-4)将C简化为设c_i和c_j是两个相邻子区域的特征值；即，其中，如果i<8,j＝i+1,如果i＝8,j＝1；简化规则为：

IFc_i＝c_j

去掉c_j

END

那么，设为中方向角所对应的空间子象限的索引向量，则

3-5)剔除不合理的设q是的长度，和是的相邻元素，其中，如果i<8,j＝i+1,如果i＝8,j＝1；则规则为：

IF q<4或存在c_i使得

停止并检索下一格点；

END

对于图3(a)实例，其是合理的。

3-6)设为中方向角所对应的空间子象限的索引向量；对步骤3-5)的和映射为单位圆的内接多边形；规则是：以或中元素的值作为单位圆的内接多边形的顶点；图3(b)的多边形映射如图3(c)和图3(d)所示。或映射能够映射的19种圆内接多边形如表1所示。

表1 19种圆内接多边形及相应形变指数

3-7)对步骤3-6)中的所有多边形构建形变指数λ，去除形变过度的候选中心，并获得确定下的涡旋的形变指数，如表1所示；λ的构建方法如下：

设圆形数据向量A＝{α₁,α₂,...,α_n}，其中，α_i≠α_j，当0°<Q(α_i,α_i+1)<180°且0°<Q(α_n,α₁)<180°时，构建形变指数λ；λ是圆形数据向量A的合成向量长度、圆形数据向量A均匀波动性和圆形数据向量A的空缺率的线性组合。

(1)合成向量长度：r_o

r_o的取值范围为[0，1]；r_o越接近1，α₁,α₂,...,α_n越聚焦于它们的均值。

(2)均匀波动性：δ_o

圆形数据向量A的有向方向差为Q(A),则Q(A)的标准差为：

同样地，表1中19类圆内接多边形的δ值中的最大值为1.73，得出标准化后的δ_o：

(3)圆形数据向量A的空缺率：ρ_o

设圆形数据向量A的维度为n,空缺率ρ_o的计算如下：

最后，将上述三种特征线性组合，得到一个综合的形变指数λ：

λ＝ar_o+bδ_o+cρ_o (9)

式(9)中，0<a<1,0<b<1,0<c<1且a+b+c＝1，a＝b＝0.25,c＝0.5，19个圆内接多边形的λ值如表1所示。

3-8)排除涡旋的伪中心点。

从风场数据中挑选涡旋样本估计λ的值：首先，每个样本的和分别由d＝5,7,9的分析域计算得到；与步骤1-3)类似，利用Bayesian决策得到最优的λ值。得到如下规则：

或

停止并检索下一格点；

END

步骤四、改进Ward聚类算法，对属于同一个涡旋系统的中心点集聚类；步骤如下：

4-1)对步骤三使用多维的窗口，即d＝5，7和9对感兴区域进行多尺度窗口的搜索候选中心点集P，其中P包含了多个涡旋系统的候选中心点集，如图4所示。利用Ward聚类算法将P中的各点聚类到各自的涡旋系统得到P_k，k＝1,2,…,K；K为类别数。

4-2)对于类别数K进行自适应选取：步骤4-1)中，每步类合并计算的最小离差平方和为D_i，计算D_i在最大离散曲率T_i处对应的类别数为最终认定的K；T_i的计算方法如下：

T_i＝1-|cosω_i| (10)

式(10)中，ω_i＝angle(p_i-1,p_i,p_i+1)是线段[p_i-1,p_i]和[p_i,p_i+1]的夹角，而K＝argmax_i{T_i}。

步骤五、利用步骤三的候选中心点集的涡度对涡旋系统分类，并精确定位全局的涡旋中心。图5分别给出了逆时针涡旋和顺时针涡旋分类和定位的详细过程。

5-1)得到了每类涡旋系统的候选中心点集P_k后，采用数学形态学中的膨胀算法对P_k进行膨胀处理，以完全覆盖涡旋中心区域，方法如下：

式(11)中，B是一个3×3的结构元素。

5-2)将涡旋类别判定为顺时针方向或逆时针方向。引进二维连续风场的垂直涡度ζ：

将式(12)转换为离散风向场的垂直涡度

式(13)中，θ_i,j是单位矢量的角，i和j是水平方向和竖直方向的索引；令l₁(或l₂)记为Z_k元素的正(或负)涡度的个数；l₁(或l₂)定义为如下：

式(14)中，#{·}集合{·}的基数；如果l₁>l₂，当前第k个涡旋判定为逆时针涡旋，否则判定为顺时针涡旋。

5-3)根据涡度值确定涡旋中心参考点。在涡旋中心区域Z_k，由公式(13)计算每个格点的垂直涡度将垂直涡度绝对值最大的检测点作为涡旋中心的参考点。

5-4)根据零风速约束精确定位涡旋中心。在到涡旋中心参考点距离小于ε的区域内，(ε＝3)，将网格划分为两个三角网格，在三角网格内通过线性插值算法计算出零风速点坐标位置，连接两个零风速点即可得到零等值线段。由此得到风速东西分量U和南北分量V在一个网格内的零等值线段，同一网格上U、V的等值线段的交点，即为风速零点。若整个局部区域存在多个风速零点，则选择距离气旋中心参考点最近的点作为气旋中心。若整个局部区域不存在风速零点，则对风速分量U、V分别进行线性插值，求解风速矢量U、V的强度极小值，将气旋中心点校正到风速矢量U、V极小值点位置。

下面以具体的实验来验证本发明提供的一种基于高精度数值风场资料的涡旋中心自动识别方法的可行性，测试样本由中国天津气象台提供，详见下文描述：

图6是利用本发明方法自动识别得到的涡旋中心定位点结果示意，该网格资料的高度为：850pha；经纬度范围为：60-150E，60N-10S；网格分辨率为：0.25°；图6中D表示逆时针涡旋中心，G表示顺时针涡旋中心，d和g代表逆时针涡旋和顺时针涡旋中心参考点。

对2014年4月含有涡旋特征的9个高精度数值风场资料进行检测验证。利用击中率POD、虚警率FAR、临界成功指数CSI和涡旋中心定准率对检验结果进行评价(见表2)。

表2检出涡旋场的临界成功指数(CSI)及涡旋中心定准率

样本	检出数	漏报数	空报数	中心校准数	击中率(％)	虚警率(％)	CSI(％)	定准率(％)
									14040808	27	1	5	24	96.4	15.6	81.8	88.9
14040908	35	1	13	30	97.2	27.1	71.4	85.7
									14041008	26	0	9	20	100	25.7	74.3	76.9
14041108	37	0	10	34	100	21.3	78.7	91.9
									14041208	41	1	14	35	97.6	25.5	73.2	85.4
14041308	41	0	13	37	100	24.1	75.9	90.2
									14041408	35	0	9	31	100	20.5	79.5	88.6
14041508	42	2	5	39	95.5	10.6	85.7	92.9
									14041608	38	0	11	34	100	22.4	77.6	89.5
总计	322	5	89	284	98.5	21.7	77.4	88.2

表2描述了9个高精度数值风场的识别情况，对其进行统计后得到涡旋样本总数为327个，成功识别样本数322个，未识别样本数5，误识别样本数89，涡旋中心校准样本数284，得出击中率、虚警率、临界成功指数和涡旋中心定准率分别为98.5％,21.7％,77.4％和88.2％。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于高精度数值风场资料的涡旋中心自动识别方法，其特征在于，包括以下步骤：

步骤一、依据850pha数值风场的矢量特征，利用衡量风向杂乱度的特征提取出包含轴对称型风场的感兴趣区域；步骤如下：

1-1)设输入的风场矩阵为W_M,N，其中，M,N分别为矩阵的行数和列数；设风场矩阵W_M,N中的两个矢量的方向夹角为Λ，构建一个统计特征即：在L×L窗口内，分别计算风场矩阵W_M,N中的网格点(i,j)和该网格点的相邻点(i，j+1)与(i-1,j)的方向夹角Λ，并累积求和；设风矢量的角度分量为θ，则两个矢量的方向夹角为Λ(θ₁,θ₂)，统计特征的计算公式如下：

${\mathrm{\&eta;}}_{m,n}^L=\underset{i=m-r}{\overset{m+r}{\&Sigma;}}\underset{j=n-r}{\overset{n+r}{\&Sigma;}}\left(\mathrm{\&Lambda;}\right({\mathrm{\&theta;}}_{i,j+1},{\mathrm{\&theta;}}_{i,j})+\mathrm{\&Lambda;}({\mathrm{\&theta;}}_{i-1,j},{\mathrm{\&theta;}}_{i,j}\left)\right)---\left(1\right)$

式(1)中，i，j，m和n分别是风场矩阵W_M,N中网格的东西方向和南北方向的索引；m＝r+2，…，M-r；n＝r+1，…，N-(r+1)，L为窗口的尺度，且为奇数；

1-2)对上述统计特征进行归一化得到归一化的统计特征对输入的给定分辨率的高精度网格风场数据，计算出统计特征的最大值，即对进行归一化得到计算公式为最后，遍历计算W_M,N，得到归一化的特征矩阵即为：

1-3)利用Bayesian决策对参数L进行优化，和对公式(1)得出的进行阈值选取，以获得输入风场矩阵的最优感兴趣区；从给定分辨率的风场数据中搜集涡旋样本V，等分为三组测试样本，包括：V_A，V_B和V_C，其中，风场数据包含不同季节、时间和海拔高度；然后，对杂乱风样本M的三个子样本：M_A，M_B和M_C和一致风样本S的三个子样本：S_A，S_B和S_C分别计算当L＝5,7,9,11和13的五类训练子特征：和当L取上述给定某一尺度值时，分别计算对于类别V_A，M_A和S_A的条件概率密度：和根据最小错误率的原则，获得的最优阈值和相应的最小错误率同样的过程应用于[V_B,M_B,S_B]和[V_C,M_C,S_C]，获得相应的最小错误率和比较不同的L取值下的最小错误率的平均值，求得对于给定分辨率的风场数据的最优L取值；基于最优L值，确定的最优阈值，该阈值包含一个与杂乱区域分界的最大值和与一致风向区域的最小值；

步骤二、利用圆形数据理论，对理想的轴对称型风场，包括：顺时针旋转场、逆时针旋转场、辐合场、辐散场和鞍形场进行特征描述和分类；步骤如下：

2-1)对于两个圆形数据变量α和β，它们之间的一种有向的方向差Q(α,β)，定义为：

定义α到β的一种走向为：

其中，α和β的取值范围各为[-180°,180°]；α到β的走向α→β为：在极坐标系下从α到β的走向；

2-2)结论1：设θ_k是以理想轴对称风场中心为中心的单位圆弧上k位置处风的角度，当0°<Q(θ_k,θ_k+ε)<180°,ε>0时，则θ_k→θ_k+ε:counterclockwise；

对于所有的0°≤k<360°，结论1都成立时，判定θ_k为第一类轴对称场，即涡旋场或辐合场或辐散场；

结论2：设θ_k是以理想轴对称风场中心为中心的单位圆弧上k位置处风的角度，当-180°<Q(θ_k,θ_k+ε)<0°,ε>0时，则θ_k→θ_k+ε:clockwise；

对于所有的0°≤k<360°，结论2都成立时，判定θ_k为第二类轴对称场，即鞍形场；

2-3)设是轴对称型风场，对于所有的0°≤k<360°，当时，和θ_k同属一类轴对称场；如果θ_k为逆时针涡旋场，当δ＝90°时，为辐合场；当δ＝180°时，为顺时针涡旋场；当δ＝270°时，为辐散场；如果θ_k为鞍形场，δ取0°～360°任意值，为鞍形场；

步骤三、相对于理想轴对称风场，实际的涡旋场是复杂的轴对称风场，依据步骤二的结论和构建变形系数，初步检索出涡旋候选中心点集；步骤如下：

3-1)对步骤一计算出的感兴趣区域的每个格点风向值进行角度量化：依次为：东风、东北风、北风、西北风、西风、西南风、南风和东南风八个方向；

3-2)以步骤一计算出的感兴趣区域为当前检索范围，依次对每个以当前格点为中心的d×d窗口下分析域的风向进行特征提取，判定当前格点是否为涡旋中心点；从中心点向X轴方向的射线开始，以45°角为步长，依次引出八条射线对分析域进行八等分；

3-3)涡旋场的特征提取；对于当前分析域的第i(i＝1,2,…,8.)个子区域，所有格点数为S_i、具有相同角度α的连通格点数为则该子区域的特征值定义为c_i，因此当前分析域的特征向量描述为C＝{c₁,c₂,…c₈}；C的提取规则为：

IF存在α使得

c_i＝α；

ELSE

停止并检索下一格点；

END

其中，ρ是阈值变量，默认值是0.4；

3-4)将C简化为设c_i和c_j是两个相邻子区域的特征值；即，其中，如果i<8,j＝i+1,如果i＝8,j＝1；简化规则为：

IFc_i＝c_j

去掉c_j

END

3-5)剔除不合理的设q是的长度，和是的相邻元素，其中，如果i<8,j＝i+1,如果i＝8,j＝1；则规则为：

IF q<4或存在c_i使得

停止并检索下一格点；

END

3-6)设为中方向角所对应的空间子象限的索引向量；对通过步骤3-5)的和映射为单位圆的内接多边形；规则是：以或中元素的值作为单位圆的内接多边形的顶点；或映射的19种圆内接多边形如表1所示；

表1

3-7)对步骤3-6)中的所有多边形构建形变指数λ，去除形变过度的候选中心，并获得确定下的涡旋的形变指数，如表1所示；λ的构建方法如下：

设圆形数据向量A＝{α₁,α₂,...,α_n}，其中，α_i≠α_j，当0°<Q(α_i,α_i+1)<180°且0°<Q(α_n,α₁)<180°时，构建形变指数λ；λ是圆形数据向量A的合成向量长度、圆形数据向量A均匀波动性和圆形数据向量A的空缺率的线性组合；

(1)合成向量长度：r_o

$r_o=\frac{1}{n}\sqrt{{\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{cos\&alpha;}}_i\right)}^2+{\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{sin\&alpha;}}_i\right)}^2}---\left(5\right)$

r_o的取值范围为[0，1]；r_o越接近1，α₁,α₂,...,α_n越聚焦于它们的均值；

(2)均匀波动性：δ_o

圆形数据向量A的有向方向差为Q(A),则Q(A)的标准差为：

$\mathrm{\&delta;}={\left\{\frac{1}{n}\mathrm{\&Sigma;}{\&lsqb;Q\left(A_i\right)-\stackrel{\&OverBar;}{Q\left(A_i\right)}\&rsqb;}^2\right\}}^{\frac{1}{2}}---\left(6\right)$

同样地，表1中所有圆内接多边形的δ值中的最大值为1.73，得出标准化后的δo：

${\mathrm{\&delta;}}_o=\frac{\mathrm{\&delta;}}{max\left\{\mathrm{\&delta;}\right\}}---\left(7\right)$

(3)空缺率：ρ_o

设圆形数据向量A的维度为n,空缺率ρ_o的计算如下：

${\mathrm{\&rho;}}_o=1-\frac{n}{8}---\left(8\right)$

最后，将上述三种特征线性组合，得到一个综合的形变指数λ：

λ＝ar_o+bδ_o+cρ_o (9)

式(9)中，0<a<1,0<b<1,0<c<1且a+b+c＝1，a＝b＝0.25,c＝0.5时，所有圆内接多边形的形变指数λ值如表1所示；

3-8)排除涡旋的伪中心点

从风场数据中挑选涡旋样本估计λ的值；首先，每个样本的和分别由d＝5,7,9的分析域计算得到；利用步骤1-3)中的Bayesian决策得到最优λ值，从而得到：

或

停止并检索下一格点；

END

步骤四、改进Ward聚类算法，对属于同一个涡旋系统的中心点集聚类；步骤如下：

4-1)对步骤三使用d＝5，7和9对感兴区域进行多尺度窗口的搜索候选中心点集P，其中P包含了多个涡旋系统的候选中心点集；利用Ward聚类算法将P中的各点聚类到各自的涡旋系统得到P_k，k＝1,2,…,K；K为类别数；

4-2)对于类别数K进行自适应选取；步骤4-1)中，每步类合并计算的最小离差平方和为D_i，计算D_i在最大离散曲率T_i处对应的类别数为最终认定的K；T_i的计算方法如下：

T_i＝1-|cosω_i| (10)

式(10)中，ω_i＝angle(p_i-1,p_i,p_i+1)是线段[p_i-1,p_i]和[p_i,p_i+1]的夹角，而K＝arg max_i{T_i}；

步骤五、利用步骤三的候选中心点集的涡度对涡旋系统分类，并精确定位全局的涡旋中心；

5-1)得到了每类涡旋系统的候选中心点集P_k后，采用数学形态学中的膨胀算法对P_k进行膨胀处理，以完全覆盖涡旋中心区域，方法如下：

Z_k＝P_k⊕B (11)

式(11)中，B是一个3×3的结构元素；

5-2)将涡旋类别判定为顺时针方向或逆时针方向；引进二维连续风场的垂直涡度ζ：

${\mathrm{\&zeta;}}_{i,j}={\&part;}_{j}v-{\&part;}_{i}u---\left(12\right)$

将式(12)转换为离散风向场的垂直涡度

$\begin{array}{c}{\widetilde{\mathrm{\&zeta;}}}_{i,j}=\frac{{\tilde{v}}_{i,j+1}-{\tilde{v}}_{i,j-1}}{2}-\frac{{\tilde{u}}_{i-1,j}-{\tilde{u}}_{i+1,j}}{2}\\ =\frac{1}{2}\&lsqb;\left({\mathrm{sin\&theta;}}_{i,j+1}+{\mathrm{cos\&theta;}}_{i+1,j}\right)-\left({\mathrm{cos\&theta;}}_{i,j-1}+{\mathrm{sin\&theta;}}_{i-1,j}\right)\&rsqb;\end{array}---\left(13\right)$

式(13)中，θ_i,j是单位矢量的角，i和j是水平方向和竖直方向的索引；令l₁(或l₂)记为Z_k元素的正或负涡度的个数；l₁或l₂定义为如下：

$l_1=\#\{(i,j)\&Element;Z_k:\widetilde{\mathrm{\&zeta;}}(i,j)>0\};l_2=\#\{(i,j)\&Element;Z_k:\widetilde{\mathrm{\&zeta;}}(i,j)<0\}---\left(14\right)$

式(14)中，#{·}集合{·}的基数；如果l₁>l₂，当前第k个涡旋判定为逆时针涡旋，否则判定为顺时针涡旋；

5-3)根据涡度值确定涡旋中心参考点；在涡旋中心区域Z_k，由公式(13)计算每个格点的垂直涡度将垂直涡度绝对值最大的检测点作为涡旋中心的参考点；

5-4)根据零风速约束精确定位涡旋中心；在到涡旋中心参考点距离小于ε的区域内，ε＝3，将网格划分为两个三角网格，在三角网格内通过线性插值算法计算出零风速点坐标位置，连接两个零风速点即可得到零等值线段；由此得到风速东西分量U和南北分量V在一个网格内的零等值线段，同一网格上U、V的等值线段的交点，即为风速零点；若整个局部区域存在多个风速零点，则选择距离气旋中心参考点最近的点作为气旋中心；若整个局部区域不存在风速零点，则对风速分量U、V分别进行线性插值，求解风速矢量U、V的强度极小值，将气旋中心点校正到风速矢量U、V极小值点位置。