一种用于气象学中自动识别并绘制冷锋的方法与流程

本发明涉及气象学领域，尤其涉及一种用于气象学中自动识别并绘制冷锋的方法。

背景技术：

冷锋是中国地区最常见的一种天气系统，它可以活动于全国各地。锋线和锋区的分析定位在天气系统中占据比较重要的地位，锋区斜压性大，有利于垂直环流的发展和能量转换。锋区可以看作风、温度、湿度等气象要素的不连续面，如能根据这些要素的特点实现冷锋的自动识别和量化描述，会有助于智能处理庞大的数值资料，挖掘客观规则，提高这类天气系统所诱发的恶劣天气的预报精准度。

随着自动识别算法在气象领域中的广泛运用，冷锋的自动识别算法不断出现。renardandclarke提出了热参数方法(tfp)的思想^[1],基于这种思想，hewson等^[2]设计了锋面自动识别的框架。这种自动识别方法得到了不少学者的认可和改进，例如：jenkner等^[3]在高分辨率条件下对此方法进行了验证，不同的是他们实验中采用了相当位温来代替湿球位温。不同于tfp思想，simmonds^[4]等利用风场漂移进行锋面的自动识别(wnd)。他们加入了锋的长度、强度等量化参数对澳洲的气象数据进行分析，验证该算法的可行性^[5]。wong^[6]利用遗传算法，根据等压线或者风场特点构造适应度函数，识别冷锋弧度片段。

发明人在实现本发明的过程中，发现现有技术中至少存在以下缺点和不足：

上述文献中所述方法大多关注一个气象要素场与冷锋的关系特征，这类方法所识别出来的锋面不尽相同，且对噪声敏感，容易造成锋面位置偏差和一些空报。国内目前还没有自动检测冷锋的方法，又由于冷锋的人工识别费时费力，会影响到预报的时效；且不能对可能出现的灾害进行及时预警，造成经济损失和人员伤亡。

[参考文献]

[1]renardr,clarkel.1965.experimentsinnumericalobjectivefrontalanalysis.mon.weatherrev.93:547–556。

[2]hewson,t.d.1998:objectivefronts.meteorologicalapplications,5(1),37-65。

[3]jenkner,j.,m.sprenger,i.schwenk,c.schwierz,s.dierer,andd.leuenberger,2010:detectionandclimatologyoffrontsinahigh-resolutionmodelreanalysisoverthealps.meteor.appl.,17,1–18。

[4]simmonds,i.,andk.keay,andj.a.t.bye,2012:identificationandclimatologyofsouthernhemispheremobilefrontsinamodernreanalysis.j.climate,25,1945–1962。

[5]hope,p.,keay,pook,m.,catto,j.,simmonds,i.,andmills,g.,2014:acomparisonofautomatedmethodsoffrontrecognitionforclimatestudies:acasestudyinsouthwestwesternaustralia.monthlyweatherreview,142(1),343-363。

[6]wong,k.y.,c.l.yip,andp.w.li.2008:automaticidentificationofweathersystemsfromnumericalweatherpredictiondatausinggeneticalgorithm.expertsyst.appl.,35,542–555。

技术实现要素：

本发明提供了一种用于气象学中自动识别并绘制冷锋的方法，本发明能自动检测出冷锋，对灾害进行及时的预警，减少经济损失和人员伤亡。

为了解决上述技术问题，本发明提出的一种用于气象学中自动识别并绘制冷锋的方法，包括以下步骤：

步骤一：基于低空风场的冷切识别；

1-1)将0°～360°的风向角进行区域划分；由大小为n×m的风向矢量阵生成用于检测冷切的三值数据阵ω3；

1-2)根据锋后为偏北风，锋前为西南风的规则，将冷切变带走向归纳为两种类型，其中，a类为自东北向西南走向，b类为自北向南走向；

检测出a类冷切变带和b类冷切变带上的点，包括：

检测出a类冷切变带上的点：对pij∈ω3，pij是三值数据阵ω3上的任意一点，圈出以pij为中心的n×n区域ωij，其中，n＜＜min{n,m}且为奇数，在三值数据阵ω3上，对ωij的左上区域ωijlu的所有点的取值求和，并记为sijlu，对ωij的右下区域ωijrd的所有点的取值求和，并记为sijrd；进而计算对a类冷切变带上的na个点的进行取值分布统计得出阈值γa，γa的取值是将的取值概率按照由大到小的顺序累加至取值分布概率和为85％时所对应的的值；定义数据阵ω4初始值为0，若则初步认为pij是该a类冷切变带上的点，并将该pij点对应的存入数据阵ω4的(i,j)位置点上；

检测出b类冷切变带上的点：对pij∈ω3，pij是三值数据阵ω3上的任意一点，圈出以pij为中心的n×m区域ωij，其中，n,m＜＜min{n,m},且为奇数，在三值数据阵ω3上，对ωij的左区域ωijl的所有点的取值求和，并记为sijl，对ωij的右区域ωijr的所有点的取值求和，并记为sijr；进而计算对b类冷切变带上的nb个点的进行取值分布统计得出阈值γb，γb的取值是将的取值概率按照由大到小的顺序累加至取值分布概率和为85％时所对应的值；若则初步认为pij是该b类冷切变带上的点，并将该pij点对应的存入数据阵ω4的(i,j)位置点上；

1-3)遍历数据阵ω4中所有点之后，将数据阵ω4中取值大于1的点置为1；

1-4)由1-2)步骤获得的检测点形成冷切变带：首先用3×3的方形结构元素对步骤1-3)获得的数据阵ω4中取值为1的点膨胀10次；然后用聚类算法将膨胀后的点聚成不同的连通区域，并计算区域面积；再使用阈值λ滤除小面积区域，λ＝1000；最后对保留区域进行骨架提取，获得待选冷切变带l；

步骤二：检测冷切变带上两侧温度差或露点差，并验证该冷切变带l是否为冷锋；

2-1)以步骤一获得的冷切变带l为基础，坐标系中，在该冷切变带l的左侧作距离为4的等距曲线l′，在该冷切变带l的右侧作距离为2的等距曲线l″，用以限定温度差或露点差的范围；运用双线性插值方法计算l′和l″上点的温度值或露点值；统计沿冷切变带l法线方向的最大温差和露点差的分布；

2-2)确定温度差阈值和露点差阈值，

根据气象预报部门提供的冷锋样本，在冷切变带l上取599个有效点统计最大温度差△tim和最大露点差△dim的分布，对599个有效点的△tim进行取值分布统计得出阈值μt，μt的取值是△tim的取值概率，按照△tim由大到小的顺序累加至取值分布概率和为75％时所对应的△tim值；对599个有效点的△dim进行取值分布统计得出阈值μd，μd的取值是△dim的取值概率，按照△dim由大到小的顺序累加至取值分布概率和为75％时所对应△dim的值；若冷切变带l上的温度差＞μt的点超过75％或露点差＞μd的点超过75％，则该冷切变带l就是冷锋；

步骤三、绘制冷锋三维示意图及其参数量化；

根据是否能够绘制出冷锋上空附近等压场的槽线，有下述情形之一：

3-1)若能够绘制出冷锋上空附近等压场的槽线，则：

3-1-1)遍历冷切变带l，得出冷锋面高空与地面夹角α的集合{αi,i＝1,...,m}，利用内嵌有可视化工具包vtk软件绘制冷锋锋面三维示意图；

3-1-2)对冷切变带l两侧的风向变化参数进行描述，变化参数包括风向最大差值、风向最小差值、风向平均差值和风向差值标准差；

3-1-3)若仅有：温度差满足冷切变带l上的温度差＞μt的点超过75％，则对温度差参数进行描述，包括温度最大差值、温度最小差值、温度平均差值和温度差值标准差；

若仅有：露点差满足冷切变带l上的露点差＞μd的点超过75％，则对露点差参数进行描述，包括露点最大差值、露点最小差值、露点平均差值和露点差值标准差；

若温度差和露点差均满足条件，则对平均差值比上标准差较大的一项参数进行描述；

3-1-4)对步骤3-1-1)冷锋锋面三维示意图中的冷锋面的倾斜角度参数进行描述，包括倾斜角度最大值、倾斜角度最小值、倾斜角度平均差值和倾斜角度标准差；

3-2)若不能绘制出冷锋上空附近等压场的槽线，则执行上述步骤3-1-2)和3-1-3)。

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

本发明提供的技术方案的有益效果是：本方法通过使用不同气象要素场上冷锋特征对冷锋进行综合识别，得到光滑的冷锋，本发明实现了冷锋的自动检测，有助于对天气灾害进行及时的预报，以减少经济损失和人员伤亡；并通过实验验证了本方法的有效性。

附图说明

图1是风向角具体区域划分方式；

图2(a)和图2(b)均为风向三值量化示例(2015年1月9日8时局部风场)，(a)风场风向分布，(b)对(a)的三值量化结果，位于图2中区域ⅰ的风向用-1表示，位于图2中区域ⅲ的风向用1表示，其他区域用0表示；

图3(a)是用于检测a类冷切变带上的点的检测范围及区域划分方案，图3(b)是用于检测b类冷切变带上的点的检测范围及区域划分方案；

图4(a)是的取值分布，横轴为图4(b)是的取值分布，横轴为纵轴为频数和频率，虚线为阈值；

图5(a)和图5(b)为两个在风向场上显示数据阵ω4中取值为1的点的运行效果实例；

图6(a)中的黑色区域为膨胀后结果，图6(b)为经阈值λ对图6(a)过滤后的结果，图6(c)为对图6(b)中的保留区域进行骨架提取的结果；

图7是关注温度差或露点差范围的说明图；

图8是双线性插值法获得格点内部点温度的示意图；

图9(a)是最大温度差分布，图9(b)是逆向累积和分布；

图10(a)是2002年06月12日08时850百帕高度场，图10(b)为02年06月04日08时850百帕高度场；其中，背景是以三角网格法得到的等压线，粗实线为自动识别的槽线；

图11(a)和图11(b)是高空槽线的地面投影和地面冷锋位置示意图；

图12是锋面三维全景示意图，时间为2015年01月06日08时，图中区域s1代表高空气压，区域s2代表锋面，区域s3为地面风场和地面温度场的叠加；

图13(a)至图13(c)是冷锋序列，右上角数字代表时间，粗实线代表击中的冷锋，虚线代表次时刻的冷锋漏识。(a)为表1中*序列，(b)为表2中^#序列，(c)为表2中*序列；

图14是本发明提供的一种用于自动识别并绘制冷锋的方法的流程图。

具体实施方式

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

本发明提出的一种用于气象学中自动识别并绘制冷锋的方法，其设计思路是：将风矢量场简化为与冷切变带相关的二值风向标记，进而转化为将无关区域标记为0的三值图像，设计点检测算法，联合膨胀、滤波以及骨架提取算法，实现冷切变带的检测、定位及切变参数提取；然后以冷切变带为基准，验证局部区域内露点差和温度差条件并给出量化参数；最后通过低压槽的识别与匹配，分析锋面的倾斜角度。本发明实现了的自动检测，有助于对天气灾害的预报，减少经济损失。同时，形成的量化描述为客观训练冷锋与极端天气的关系模型奠定了一定基础。

如图14所示，该方法主要包括：基于低空风场的冷切识别，冷切变带上两侧温度差或露点差检测及验证，冷切变带上空附近等压场上槽线检测和冷锋的三维示意及量化参数；具体内容如下：

步骤一：基于低空风场的冷切识别；

1-1)将0°～360°的风向角进行区域划分；具体区域划分方式如图1所示，区域ⅰ意于包含冷锋锋前西南风的风向范围，α1和α2表示对西南风的放松程度，本发明中，令α1＝α2＝10°；区域ⅲ意于包含冷锋锋后的偏北风的风向范围，α3表示对西北风的放松程度，α4表示对正北风的放松程度，本发明中，令α3＝20°，α4＝45°。区域ⅱ和区域ⅳ所示的风向范围与冷切无关。

1-2)由大小为n×m的风向矢量阵生成用于检测冷切的三值数据阵ω3；即qij是风向矢量阵上的任意一点，θij是qij的风向角，根据式(1)算出三值标记aij，且将三值标记aij存于三值数据阵ω3中，

图2(a)和图2(b)给出了经式(1)处理的三值量化的示意图，其中，图2(a)为风场风向分布，图2(b)是对图2(a)的三值量化结果。

1-2)根据锋后为偏北风，锋前为西南风的规则，将冷切变带走向归纳为两种类型，其中，a类为自东北向西南走向，b类为自北向南走向；

检测出a类冷切变带和b类冷切变带上的点，包括：

检测出a类冷切变带上的点：对pij∈ω3，pij是三值数据阵ω3上的任意一点，圈出以pij为中心的n×n区域ωij，其中，n＜＜min{n,m}且为奇数，图3(a)为用于检测a类冷切变带上的点的检测范围及区域划分方案。在三值数据阵ω3上，对图3(a)中左上区域ωijlu的所有点的取值求和，并记为sijlu，对图3(a)中的右下区域ωijrd的所有点的取值求和，并记为sijrd；进而计算对a类冷切变带上的na个点的进行取值分布统计得出阈值γa，γa的取值是将的取值概率，按照由大到小的顺序累加至取值分布概率和为85％时所对应的的值，统计结果如图4(a)所示，可知γa＝34；定义数据阵ω4初始值为0，若则初步认为pij是该a类冷切变带上的点，并将该pij点对应的存入数据阵ω4的(i,j)位置点上；

检测出b类冷切变带上的点：对pij∈ω3，pij是三值数据阵ω3上的任意一点，圈出以pij为中心的n×m区域ωij，其中，n,m＜＜min{n,m},且为奇数，图3(b)为用于检测b类冷切变带上的点的检测范围及区域划分方案。在三值数据阵ω3上，对图3(b)中左区域ωijl的所有点的取值求和，并记为sijl，对图3(b)中右区域ωijr的所有点的取值求和，并记为sijr；进而计算对b类冷切变带上的nb个点的进行取值分布统计得出阈值γb，γb的取值是将的取值概率，按照由大到小的顺序累加至取值分布概率和为85％时所对应的值,统计结果如图4(b)所示，可知γb＝22；若则初步认为pij是该b类冷切变带上的点，并将该pij点对应的存入数据阵ω4的(i,j)位置点上；

1-3)遍历数据阵ω4中所有点之后，将数据阵ω4中取值大于1的点置为1；图5(a)和图5(b)为两个在风向场上显示数据阵ω4中取值为1的点的运行效果实例。

1-4)由1-3)步骤获得的检测点形成冷切变带：首先用3×3的方形结构元素对步骤1-3)获得的数据阵ω4中取值为1的点膨胀10次；然后用聚类算法将膨胀后的点聚成不同的连通区域，并计算区域面积；再使用阈值λ滤除小面积区域，λ＝1000；最后对保留区域进行骨架提取，获得待选冷切变带l；图6(a)中的黑色区域为膨胀后结果,图6(b)为经阈值λ对图6(a)过滤后的结果，图6(c)为对图6(b)中的保留区域进行骨架提取的结果。

步骤二：检测冷切变带上两侧温度差或露点差，并验证该冷切变带l是否为冷锋；

2-1)以步骤一获得的冷切变带l为基础，在该冷切变带l的左侧作距离为4(即坐标系中的距离为4，而实际距离为100km)的等距曲线l′，在该冷切变带l的右侧作距离为2(坐标系中距离为2，而实际距离为50km)的等距曲线l″，用以限定温度差或露点差的范围；

图7中曲线l是已检出的冷切变带，设b是冷切变带l上的任意一点，a、c是其邻点，b点的切线方向用线段ac的斜率来估计；由属于l′的点d和属于l″的点e得到线段de，de垂直于b点的切线方向，面对分辨率为0.25×0.25的气象格点数据，d、b间的距离为rdb＝4，e、b间的距离为reb＝2。为计算等距曲线l′和等距曲线l″上各点的坐标，首先计算d点和e点的坐标，然后算出d点和e点相对于b点的偏移量。鉴于检测出的切变线呈现单调性，d点相对于b点的偏移量即为等距曲线l′对冷切变带l的偏移量，即冷切变带l的点加上此偏移量就能得到等距曲线l′相应点的坐标；e点相对于b点的偏移量即为等距曲线l″对冷切变带l的偏移量，即冷切变带l的点加上此偏移量就能得到等距曲线l″相应点的坐标。

设：a～e的坐标依次为(xa,ya)、(xb,yb)、(xc,yc)、(xd,yd)和(xe,ye)

则：ac两点间线段的斜率：

b点的法线方向：

de所在直线的方程：y＝k′x+b(4)

该直线上的任意一点到b点的距离：

对式(4)和式(5)进行联立求解，得到d点坐标(xd,yd)和e点坐标(xe,ye)，进而得出图7中的(△x1,△y1)和(△x2,△y2)。对于等距曲线l′上的任意一点，其坐标(xd′，yd′)可由冷切变带l上对应的(xb′，yb′)根据下式计算得出。同理，对于等距曲线l″上的任意一点，其坐标(xe′，ye′)可由冷切变带l上对应的(xb′，yb′)也可根据下式计算得出。

2-2)运用双线性插值方法计算l′和l″上点的温度值或露点值；统计沿冷切变带l法线方向的最大温差和露点差的分布；下文以温度为例，给出具体计算方法。

设d为l′任意一点且未落于格点上，并设与d点相邻的四个格点分别为m～q，如图8所示，坐标依次为(xm,ym)、(xn,yn)、(xp,yp)和(xq,yq)，它们的温度值依次为tm、tn、tp和tq，利用线性插值方法不难算出d0点温度值td0和d1点温度值td1：

利用td0和td1再次线性插值得到d点温度值td：

计算其法线方向上的温度值，进而计算出最大温度差值。

2-3)确定温度差阈值和露点差阈值，

根据气象预报部门提供的冷锋样本，在冷切变带l上取599个有效点统计最大温度差△tim和最大露点差△dim的分布，对599个有效点的△tim进行取值分布统计得出阈值μt，μt的取值是△tim的取值概率，按照△tim由大到小的顺序累加至取值分布概率和为75％时所对应的△tim值；图9(a)给出了599个点的最大温度差△tim分布，由最大温度差值开始，沿图9(a)的横轴逆向累积求和，获得一条单调上升累计和曲线，如图9(b)所示，据此，选择曲线开始变缓且和数达到75％的点作为阈值μt(＝4℃)，同样方法，得到冷切变带前后露点差阈值μd(＝4℃)。若冷切变带l上的温度差＞μt的点超过75％或露点差＞μd的点超过75％，则该冷切变带l就是冷锋；

步骤三、绘制冷锋三维示意图及其参数量化；

根据是否能够绘制出冷锋上空附近等压场的槽线，有下述情形之一：

3-1)若能够绘制出冷锋上空附近等压场的槽线，则：

3-1-1)遍历冷切变带l，得出冷锋面高空与地面夹角α的集合{αi,i＝1,...,m}，α的计算方法如下所示：

图10(a)是2002年06月12日08时850百帕高度场，图10(b)为02年06月04日08时850百帕高度场；其中，背景是以三角网格法得到的等压线，粗实线为自动识别的槽线。图11(a)中，曲线s′是高空槽线s在地面上的投影。设b点是冷切变带l上的任意一点，a′点是b点在曲线s′上的投影点(若无投影点则舍弃此点)，图11(b)为高空槽线s上点a和地面冷切变带上b点沿a-b形成的垂直剖面示意图。

设700hpa等压面到地面距离为h，a′点和b点坐标分别为(xa′,ya′)和(xb,yb)，则夹角α：

其中，

利用内嵌有可视化工具包vtk(visualizationtoolkit)的visualstudio2010软件绘制冷锋锋面三维示意图，如图12所示。

3-1-2)对冷切变带l两侧的风向变化参数{△θi,i＝1,...,n}进行描述，变化参数包括风向最大差值θm、风向最小差值θm、风向平均差值和风向差值标准差σθ，记作

3-1-3)若只有：温度差满足冷切变带l上的温度差＞μt的点超过75％这一条件，则对温度差参数{△tim,i＝1,...,n}进行描述，包括温度最大差值△tmm、温度最小差值△tmm、温度平均差值和温度差值标准差σ△tm，记作

若只有露点差满足：冷切变带l上的露点差＞μd的点超过75％这一条件，则对露点差参数{△dim,i＝1,...,n}进行描述，包括露点最大差值△dmm、露点最小差值△dmm、露点平均差值和露点差值标准差σ△dm，记作

若温度差和露点差均满足条件，则对平均差值比上标准差较大的一项参数进行描述；

3-1-4)对步骤3-1-1)冷锋锋面三维示意图中的冷锋面的倾斜角度参数{αi,i＝1,...,m}进行描述，包括倾斜角度最大值αm、倾斜角度最小值αm、倾斜角度平均差值和倾斜角度标准差σα，记作根据均值将冷锋记为前倾锋(槽线在冷锋前面)，后倾锋(槽线在冷锋后面)；

槽线的绘制可以根据专利申请号为2015103197143的专利文献中披露的技术方案来实现：首先采用断面极值法进行槽线特征点分析，然后根据槽线的走向，将梯度方向为-180°、-135°和-90°的气压格点数据作为槽线候选特征点；再对槽线候选特征点进行纠偏、过滤掉不满足条件的槽线候选点；通过光滑中轴算法绘制出槽线。

3-2)若不能绘制出冷锋上空附近等压场的槽线，则执行上述步骤3-1-2)和3-1-3)。

下面以具体的测试来验证本发明实施例提供的一种冷锋系统的识别方法的可行性，详见以下描述：

数据集1：20150101-20150120的早8点和晚8点共20组数据，该系列数据参与了算法参数的调试。

数据集2：20161220-20161229的早8点和晚8点共20组数据，该系列数据未参与算法参数的调试。

表1和表2分别列出了本发明方法对数据集1和数据集2的识别结果。√表示是，×表示否，“空”代表此步骤不执行，“击中”表示存在冷锋并且正确识别，“空报”表示不存在冷锋但是有识别结果，“漏报”表示存在冷锋但是没有给出识别结果，“正确”表示不存在冷锋并且没有识别结果。其中，数据集1检出3个冷锋过程，数据集2检测出2个冷锋过程(122008时刻中一个击中的冷锋为该时刻之前系统的延续，下文的结果统计忽略此次冷锋)。图13(a)至图13(c)展示了其中三个时长超过2天的冷锋及其移动情况。表3是图13(a)至图13(c)中冷锋的量化参数，空代表无高空槽线。

表1数据集1实际检出识别结果统计(2015年)

表2数据集2实际检出识别结果统计(2016年)

表3量化描述

数据集1中实际持续存在3个冷锋系统，本发明方法均持续击中，未出现漏识和空报；数据集2中实际持续存在2个冷锋系统，本发明方法各漏识1个时次，共空报1个时次。

在对数据集1和数据集2运行本发明算法的同时，对各步耗时进行了统计，表4列出了它们的平均运行时间。

表4算法平均运行时间

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。