CN104951624A - 一种计算机气象软件中基于风场数据的槽线自动绘制方法 - Google Patents

Abstract

本发明公开了一种计算机气象软件中基于风场数据的槽线自动绘制方法，实现天气图中气象槽线的自动识别，包括以下步骤：长波槽生成：提取风场网格数据的逆向环流中心，查找数据中的候选长波槽点，以逆向环流中心为起点进行槽点追踪和连接，提取长波槽。短波槽生成：对风场网格数据进行涡度、角度偏转筛选，提取候选短波槽点，在结合候选长波槽点进行筛选后，再对余下候选短波槽点进行最小生成树聚类，对各个生成类中的候选槽点集分别进行曲线拟合，提取短波槽。

Description

一种计算机气象软件中基于风场数据的槽线自动绘制方法

技术领域

本发明涉及气象槽线自动分析，属于气象分析中的自动化分析领域，具体地说是一种计算机气象软件中基于风场数据的槽线自动绘制方法。

背景技术

天气图是填有各地同一时间气象观测记录的特种地图，它描述了某一瞬间某一区域的天气情况。在日常的气象业务工作中，需要绘制地面天气图(简称地面图)和等压面图(简称高空图)两类。其中地面天气图填有地面各种气象要素和提拿起现象，如气温、湿度、风向、风速、海平面气压、雨、雪、雾和云高云状等。等压面图是用来表示等压面起伏形势的图，常用的等压面绝对形势图有：850hpa，700hpa，500hpa等压面图。等压面图分析的主要内容包括等高线的分析，等温线的分析，湿度场的分析，槽线和切变线的分析，温度平流和湿度平流。其中已经实现自动分析的功能包括等值线的自动绘制。而槽线作为高空图中很重要的一项还没有实现自动分析。

一般来说槽线与脊线相对应，它是等值线曲率最大点按一定规则的连线，在地形数据分析过程中称为山脊线，山谷线。在气象领域，槽线定义为低压槽区内等高线曲率最大点的连线，是气旋中心到低压槽内气压最低点的集合所组成的曲线。在天气分析中，槽线的分析有着非常重要的作用，它一定程度上可以反映出天气的变化趋势，可以为预报人员提供重要的参考，例如，在水平方向上槽线的前方多为西南暖湿气流，后方盛行干燥的西北气流，而暖湿气流通常可带来降雨，干燥气流一般使得天气晴朗；在垂直方向上，槽线前方有气流的上升运动，水汽充沛，可产生降水，后方则一般为气流下沉伴随天气转晴。综上所述，槽线是预报中成云致雨的重要依据。

准确判断槽线是天气系统分析和天气图绘制的重要内容。美国AWIPS系统与国家气象局MICAPS系统均已实现天气图的自动填绘和等值线的自动绘制等功能，但都没有实现槽线的自动分析。目前，气象预报中的槽线分析基本上采用人工分析的方法实现。人工分析槽线存在着明显的缺陷，一方面分析工作量较大，效率较低；另一方面对于槽线的分析，预报员主观上分析经验不同，客观上侧重依据差异，都对槽线的分析产生影响，不利于改进提高及应用推广。因此，建立一个稳定的、准确的系统来实现槽线的自动绘制，对于实现气象自动分析、提高工作效率具有重要的理论研究意义和实际应用价值。

发明内容

发明目的：本发明所要解决的技术问题是针对目前槽线自动分析不准确、不完善，提出了一种计算机气象软件中基于风场数据的槽线自动绘制方法。

技术方案：本发明公开一种基于风场网格数据的槽线自动识别方法，包括以下步骤：

步骤1，长波槽生成：对风场网格数据进行特征点提取，判别逆向环流中心，同时提取候选长波槽点，对于每一个逆向环流中心，将其作为起点进行槽点追踪，输出单条长波槽，并将已使用槽点标注；

步骤2，短波槽生成：标记步骤1中长波槽区域，对风场网格数据进行涡度、角度偏转筛选得到候选短波槽点，排除长波槽区域内的候选短波槽点后，对剩下的候选短波槽点的集合进行最小生成树聚类，在各个生成类的候选短波槽点集合中进行曲线拟合作为短波槽分析结果；

步骤3，计算机气象软件根据步骤2的结果，自动绘制风场的槽线并显示。

其中，步骤1包括以下步骤：

步骤1-1，提取逆向环流中心：在风场网格数据中查找局部风速最小点得到特征点，计算特征点处雅可比Jacobian矩阵并判断矢量拓扑结构，得到逆向环流中心；

步骤1-2，提取候选长波槽点：遍历风场网格数据，计算网格点风矢方向，依据水平、垂直两个方向上相邻网格风矢方向特点判断候选槽点，得到候选槽点集合；

步骤1-3，追踪槽线：通过步骤1-1、步骤1-2分别得到逆向环流中心和候选槽点集合，将用于南向追踪的候选槽点记为WP，WP构成的候选槽点集合为S_W，将用于东向追踪的候选槽点记为EP，EP构成的候选槽点集合为S_E，由长波槽起点开始分两个方向分别追踪两类候选槽点，分别得到槽线分析结果为点列PList_W和PList_E，即单条长波槽。

步骤1-1包括以下步骤：

步骤1-1-1，查找特征点：遍历风场网格数据，在每一个3*3网格中心查找局部风速最小点，选取对称位置上风向相反的点作为特征点；

步骤1-1-2，筛选逆向环流中心：构建Jacobian矩阵，求矩阵特征值，通过特征值判断特征点类型，筛选出所需逆向环流中心。

步骤1-1-1包括以下步骤：

步骤1-1-1-1，查找局部风速最小点：从风场网格数据左上角开始选取一个网格数据3*3点阵，寻找局部风速最小点，该点阵形式如下表所示：

V11	V12	V13
			V21	V22	V23
V31	V32	V33

每个方格代表风场中的一个网格点，V₁₁,V₁₂,...V₃₃代表在各自网格点的风速大小，由下式计算：

$V_{ii,jj}=\sqrt{u_{ii,jj}^2+v_{ii,jj}^2},$

其中，ii、jj分别为网格点行、列坐标索引，ii、jj取值范围1～3，u_ij,v_ij分别为网格点坐标(ii,jj)处的风速水平分量与垂直分量，在该点阵中，如果中心点V₂₂的值比它周围的8个点处的值都小，则判别该点为局部风速最小点，并将其作为一个特征点候选点；

步骤1-1-1-2，结合风向判断：对于满足步骤1-1-1-1的特征点候选点，判断关于该特征点候选点成中心对称的网格点上的风速方向是否相反，即在该3*3矢量点阵中，是否满足以下条件：

u₁₁×u₃₃＜0且v₁₁×v₃₃＜0，

或

u₁₃×u₃₁＜0且v₁₃×v₃₁＜0，

如满足，该候选点即为特征点，否则判定不是特征点；

步骤1-1-1-3，按行列顺序坐标增1对风场网格数据所有3*3点阵，重复步骤1-1-1-1～1-1-1-2直到最后一个网格结束；

步骤1-1-2包括以下步骤：

步骤1-1-2-1，求解Jacobian矩阵J：

$J=\left(\begin{array}{cc}\frac{\∂u}{\∂x}& \frac{\∂u}{\∂y}\\ \frac{\∂v}{\∂x}& \frac{\∂v}{\∂y}\end{array}\right),$

其中u,v分别为任意一点处风速水平分量与垂直分量，分别表示该点风速水平分量u在水平与垂直方向上的偏导；分别表示该点风速垂直分量v在水平与垂直方向上的偏导。

对于矩阵中非边界数据点，即下标大于1，行标小于矩阵行数，列表小于矩阵列数的网格点，坐标索引为(i,j)的网格点处偏导用下式表示：

$\frac{\∂u_{i,j}}{\∂x}=\frac{1}{2}(u_{i,j+1}-u_{i,j-1})\frac{\∂u_{i,j}}{\∂y}=\frac{1}{2}(u_{i-1,j}-u_{i+1,j}),$

$\frac{\∂v_{i,j}}{\∂x}=\frac{1}{2}(v_{i,j+1}-v_{i,j-1})\frac{\∂v_{i,j}}{\∂y}=\frac{1}{2}(v_{i-1,j}-v_{i+1,j}),$

首行首列点处的偏导求解采取下式表示：

$\frac{\∂u_{i,1}}{\∂x}=u_{i,2}-u_{i,1}\frac{\∂u_{1,j}}{\∂y}=u_{1,j}-u_{2,j},$

$\frac{\∂v_{i,1}}{\∂x}=v_{i,2}-v_{i,1}\frac{\∂v_{1,j}}{\∂y}=v_{1,j}-v_{2,j},$

末行Row、末列Col点处的偏导求解采用下式表示：

$\frac{\∂u_{i,Col}}{\∂x}=u_{i,Col}-u_{i,Col-1}\frac{\∂u_{Row,j}}{\∂y}=u_{Row-1,j}-u_{Row,j},$

$\frac{\∂v_{i,Row}}{\∂x}=v_{i,Row}-v_{i,Row-1}\frac{\∂v_{Col,j}}{\∂y}=v_{Col-1,j}-v_{Col,j},$

步骤1-1-2-2，求解特征值：采用下式求解坐标索引为(i,j)的网格点处Jacobian矩阵J_i,j：

$J_{i,j}=\left(\begin{array}{cc}\frac{\∂u_{i,j}}{\∂x}& \frac{\∂u_{i,j}}{\∂y}\\ \frac{\∂v_{i,j}}{\∂x}& \frac{\∂v_{i,j}}{\∂y}\end{array}\right),$

求解行列式：

$\left|\begin{array}{cc}\frac{\∂u_{i,j}}{\∂x}-\mathrm{\λ}& \frac{\∂u_{i,j}}{\∂y}\\ \frac{\∂v_{i,j}}{\∂x}& \frac{\∂v_{i,j}}{\∂y}-\mathrm{\λ}\end{array}\right|=0,$

得到该点的Jacobian矩阵的特征值，即λ的两个解，记为λ₁,λ₂；

步骤1-1-2-3，判断特征点类型：根据矢量拓扑结构特性，提取其中吸引旋转点和旋转点两种拓扑结构即为逆向环流中心：

通过步骤1-1-2-1计算得到的特征值λ₁,λ₂，其中λ₁对应的实部与虚部分别为R₁,I₁，λ₂对应实部和虚部分别为R₂,I₂，当R₁＜0,R₂＜0；I₁≠0,I₂≠0时，判定特征点为吸引旋转点；当R₁＝0,R₂＝0；I₁≠0,I₂≠0时，判定特征点为旋转点；

步骤1-1-2-4，对于每一个特征点，重复步骤1-1-2-1～1-1-2-3直到最后一个特征点，得到逆向环流中心点集合。

步骤1-2包括以下步骤：

步骤1-2-1，计算风矢方向：在风场网格数据中，建立弧度坐标系，北风弧度为9东°，东风的弧度为180°，南风弧度为270°。通过风矢水平、垂直两个方向上的分量求得每一个网格点处风矢的方向。

风场网格数据的网格尺寸为Row*Col，P_i,j表示坐标索引为(i,j)的网格点，P_i,j.X为P_i,j点处X方向上的几何坐标，P_i,j.Y为P_i,j点处Y方向上的几何坐标，P_i,j处风矢方向为α_i,j，分别计算一个由四个相邻网格点组成的网格中各点处弧度，即计算P_i,j、P_i,j+1、P_i+1,j、P_i+1,j+1处的方向为α_i,j、α_i,j+1、α_i+1,j、α_i+1,j+1；

步骤1-2-2，判断α_i,j与α_i,j+1角度范围：若α_i,j＜90°且270°＜α_i,j+1＜360°，则P_i,j、P_i,j+1两点间存在用于南向追踪的候选槽点WP，候选槽点WP在X方向上、Y方向上的几何坐标WP.X和WP.Y分别为：

$WP.X=\frac{180-{\mathrm{\α}}_{i,j}}{{\mathrm{\α}}_{i,j+1}-{\mathrm{\α}}_{i,j}}*(P_{i,j+1}.X-P_{i,j}.X)+P_{i,j}.X,$

WP.Y＝P_i,j.Y；

步骤1-2-3，判断α_i,j与α_i+1,j角度范围：若180°＜α_i,j＜270°且270°＜α_i+1,j＜360°，那么P_i,j、P_i+1,j两点间存在用于东向追踪的候选槽点EP，候选槽点EP在X方向上、Y方向上的几何坐标EP.X和EP.Y分别为：

EP.X＝P_i,j.X，

$EP.Y=\frac{180-{\mathrm{\α}}_{i,j}}{{\mathrm{\α}}_{i+1,j}-{\mathrm{\α}}_{i,j}}*(R_{i+1.j}.Y-P_{i,j}.Y)+P_{i,j}.Y;$

步骤1-2-4，遍历风场网格中所有网格点数据，重复步骤1-2-1～1-2-3筛选候选长波槽点直到最后一个网格，得到候选槽点集合。

步骤1-3包括以下步骤：

步骤1-3-1，选择追踪起点：从逆向环流中心集合中按序任意选取一个点作为长波槽追踪起点StartP，当追踪方向Direction＝270°时，设候选槽点集合S＝S_W，点列PList＝PList_W；当追踪方向Direction＝180°，设候选槽点集合S＝S_E，点列PList＝PList_E；

步骤1-3-2，筛选候选长波槽点：通过距离、方向筛选选取符合条件的候选槽点：

步骤1-3-2-1，距离筛选：遍历候选槽点集合S，通过距离阈值ThredL对其进行筛选，设水平相邻网格点间距为disAdjacent，网格点间距是经纬度上网格点风矢对应到绘制几何坐标系中两个水平相邻网格点间的几何距离，由绘制风矢时选取的距离决定。ThredL应为3*disAdjacent～5*disAdjacent之间。点P_k为S集合中第k点，点P_k与追踪起点StartP距离为L_k，若L_k＜ThredL，将P_k加入筛选后候选槽点集合S1中，否则舍弃，若不存在满足L_k＜ThredL的点，返回执行步骤1-3-1；

步骤1-3-2-2，方向筛选：遍历候选槽点集合S1，点P_m为候选槽点集合S1中第m点，点P_m与追踪起点StartP构成的直线与追踪方向Direction夹角为γ_m，

候选槽点与追踪方向构成的弧度满足以下弧度偏转的条件：

若槽线已探索延伸至点P₁，追踪方向为AB，下一个可能在槽线上的关键点应在以AB为轴一定大小的张角范围内。那么左侧存在一个可偏转的弧度，其阈值设定为α_left，右侧存在一个可偏转弧度阈值，其阈值设定为α_right，设置α_left＝30°，α_right＝30°。槽线的追踪设置在上述规定角度范围内；将满足上述条件且γ_m＜Thredγ的P_m加入筛选后的候选槽点集合S2，若不存在满足上述条件的点，返回步骤1-3-1；

步骤1-3-3，通过评估方程确定下一候选点：设置评估方程Evaluation(L,α)对每一个符合条件的候选槽点进行评价：

$Evaluation(L,{\mathrm{\α}}_{field})=\frac{WeightsA}{L^2}+\frac{WeightsB}{|{\mathrm{\α}}_{field}-{\mathrm{\α}}_{fixed}|},$

其中，L为候选槽点与已延伸槽线上最后一点的距离，α_field为候选槽点与最后一点构成线段与追踪方向的夹角，α_fixed为使得槽线追踪向着需要的方向偏转的固定偏转量，WeightsA和WeightsB分别为距离判断和弧度偏转判断权值，这两个权值的设定是对距离和弧度偏转两个要素的重要性进行权衡，取值范围均为0～1。遍历候选槽点集合S2，查找使得Evaluation(L,α_field)最大的点P_Next，将点P_Next加入点列PList；

步骤1-3-4，重确定追踪起点、追踪方向：求取追踪起点StartP与下一个点P_Next在坐标系中角度δ，设置追踪方向Direction＝δ+α_fixed，将追踪起点StartP设置为P_Next；

步骤1-3-5，按序遍历逆向环流中心集合，重复步骤1-3-2、1-3-3、1-3-4直到最后一个逆向环流中心点。

步骤2包括以下步骤：

步骤2-1，标记长波槽区域：由步骤1-3得到组成长波槽的点列，对这些点周围一定范围进行标记，该范围距离应该保证在4*disAdjacent～6*disAdjacent之间，disAdjacent由步骤1-3-2-1得到。得到标记的区域M；

步骤2-2，筛选候选短波槽点：包括如下步骤：

步骤2-2-1，筛选偏转角度：对风场网格数据中的每一个网格点，分别在水平和垂直两个方向上求取相邻网格点风矢的方向偏转大小，包括如下步骤：

步骤2-2-1-1，求取角度偏转量：Angle_horizon_s,t、Angle_vertical_s,t分别为网格点(s,t)处水平与垂直方向上的角度偏转，则：

Angle_horizon_s,t＝α_s,t-α_s,t+1，

Angle_vertical_s,t＝α_s,t-α_s+1,t，

若Angle_horizon_s,t＜0，设置Angle_horizon_s,t＝-Angle_horizon_s,t；

若Angle_vertical_s,t＜0，设置Angle_vertical_s,t＝-Angle_vertical_s,t；

若Angle_vertical_s,t＞180°，设置Angle_vertical_s,t＝360°-Angle_vertical_s,t；

若Angle_horizon_s,t＞180°，设置Angle_horizon_s,t＝360°-Angle_horizon_s,t；

步骤2-2-1-2，完成步骤2-2-1-1角度偏转计算后，计算偏转量的筛选阈值AngleThred：

$\begin{array}{c}AngleThred=\mathrm{\β}(\frac{1}{Col*(Row-1)}\underset{s,t}{\Σ}Angle\_{\mathrm{horizon}}_{s,t}\\ +\frac{1}{(Col-1)*Row}\underset{s,t}{\Σ}Angle\_{\mathrm{vertical}}_{s,t})\end{array},$

其中β为可调整权值，取值范围为0～5；

步骤2-2-1-3，遍历所有网格点，若Angle_vertical_s,t＞AngleThred，则将网格中第s行、第t列点与第s行、第t+1列点之间中点标记为候选槽点，加入候选短波槽点集合O₁；若Angle_horizon_s,t＞AngleThred，则将网格中第s行、第t列点与第s+1行、第t列点之间中点标记为候选槽点，加入候选短波槽点集合O₁；

步骤2-2-2，筛选涡度：遍历风场网格中网格点数据求得各点处涡度值，计算涡度筛选阈值进行筛选提取候选短波槽点：

步骤2-2-2-1，在风场网格数据中，网格点(i,j)处涡度大小Vorticity_i,j为：

Vorticity_i,j＝|u_x-v_y|

其中u_x是水平分量u在水平方向上的偏导，v_y是垂直分量v在垂直方向上的偏导，同步骤1-1-2-1用差分代替上式中的偏微分计算涡度大小；

步骤2-2-2-2，计算涡度筛选阈值VorticityThred：

$VorticityThred=\mathrm{\γ}*\left(\frac{1}{Row*Col}\underset{i,j}{\Σ}{\mathrm{Vorticity}}_{i,j}\right),$

其中，γ是可调整的涡度筛选权重，取值范围为0～5；

步骤2-2-2-3，标记候选短波槽点：若Vorticity_i,j＞VorticityThred，则将风场网格中第i行、第j列点标记为候选短波槽点，将其加入候选短波槽点集合O₁；

步骤2-3，排除部分候选短波槽点：结合步骤2-1已标记的部分区域M，筛选后的候选槽点集合为O₂：

O₂＝O₁-(O₁∩M)；

步骤2-4，候选槽点最小生成树聚类：对候选槽点集合O₂中的点依据点集合坐标位置建立加权无向连通图，进行最小生成树聚类，生成多个分类：

步骤2-5，曲线拟合：对生成的分类设定拟合多项式，通过最小二乘法确定参数，实现多项式曲线拟合，提取短波槽分析结果。

步骤2-4包括以下步骤：

步骤2-4-1，建立加权无向连通图G＝(Q,E)：将候选槽点集合O₂中的槽点按需排列编号为q₁,q₂,...,q_n，n为候选槽点集合O₂中点的个数，顶点集合为Q＝{q₁,q₂,...,q_n}，边集合为E，E＝{＜q₁,q₂＞,＜q₁,q₃＞,...,＜q_n-1,q_n＞}，其中＜q_i,q_j＞为顶点q_i与q_j的距离；

步骤2-4-2，初始化顶点集合：初始化顶点集合Q_new＝{q_z}，其中q_z为Q中任意一点，z取值范围为自然数，初始边集合E_new＝{}为空；

步骤2-4-3，构建最小生成树：从初始化顶点集合Q_new＝{q_z}中的起点开始，在边集合E中查找与顶点集合Q_new中的点构成的边中权值最小的边，如果这条边上的点不在顶点集合Q_new中，则将该点加入到顶点集合Q_new中：

步骤2-4-3-1，在边集合E中选取权值最小的边＜q_l,q_d＞，其中q_l为顶点集合Q_new中的元素，q_d不在顶点集合Q_new中，如果存在有多条边权值均为最小值，则任意选取其中之一；

步骤2-4-3-2，将q_d加入到顶点集合Q_new中，将＜q_l,q_d＞边加入到边集合E_new中；

步骤2-4-3-3，重复步骤2-4-3-1～2-4-3-2直到Q_new＝Q，输出聚类结果，使用集合Q_new、E_new描述所得最小生成树；

步骤2-4-4，分割子树：对步骤2-4-3-3得到的最小生成树边集合E_new，设置筛选阈值ClustThred，，ClustThred应在3*disAdjacent～5*disAdjacent范围内。遍历边集合E_new所有元素，若其权值大于ClustThred，则去掉该条边，将最小生成树分为两个以上的分树，即可得到候选短波槽点集合通过最小生成树方法划分的子类。

步骤2-5包括以下步骤：

步骤2-5-1，设定多项式：依据多项式f(h)＝a₁h⁴+a₂h³+a₃h²+a₄h+a₅进行拟合，其中h为自变量，f(h)为因变量，a₁,a₂,a₃,a₄,a₅为待确定参数，对一个总共g个点的分类点集合Sub＝{p₁,p₂,...,p_g}，g取值范围为自然数，每一个点由水平方向上和垂直方向上两个坐标决定，分别设为p_r.X和p_r.Y，r＝1,2,...,g，将这些点水平方向上的坐标作为自变量，垂直方向上的坐标作为因变量；

步骤2-5-2，采用最小二乘法求取参数：依据偏差平方和最小的原则求取待确定参数，

距离平方和即

分别对等式右侧a₁,a₂,a₃,a₄,a₅求偏导得：

令上式都为零得到范德蒙德矩阵：

$\left[\begin{array}{ccccc}g& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^7\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^7& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^8\end{array}\right]*\left[\begin{array}{c}{a}_5\\ a_4\\ a_3\\ a_2\\ a_1\end{array}\right]=\left[\begin{array}{c}\underset{r=1}{\overset{g}{\Σ}}{p}_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4*p_r.Y\end{array}\right],$

将该范德蒙德矩阵简化后得到：

$\left[\begin{array}{ccccc}1& p_1.X& p_1.X^2& p_1.X^3& p_1.X^4\\ 1& p_2.X& p_2.X^2& p_2.X^3& p_2.X^4\\ 1& p_3.X& p_3.X^2& p_3.X^3& p_3.X^4\\ .& .& .& .& .\\ .& .& .& .& .\\ 1& p_g.X& p_g.X^2& p_g.X^3& p_g.X^4\end{array}\right]*\left[\begin{array}{c}{a}_5\\ a_4\\ a_3\\ a_2\\ a_1\end{array}\right]=\left[\begin{array}{c}{p}_1.Y\\ p_2.Y\\ p_3.Y\\ .\\ .\\ p_g.Y\end{array}\right],$

求得a₁,a₂,a₃,a₄,a₅的值，即得多项式拟合曲线，将其作为各分类点集合短波槽分析结果。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和或其他方面的优点将会变得更加清楚。

图1为本发明的总体流程示意图。

图2为长波槽提取过程示意图。

图3为短波槽提取过程示意图。

图4a和图4b为风矢绘制结果。

图5为提取逆向环流中心示意图。

图6为提取长波候选槽点示意图。

图7为对长波槽进行追踪得到追踪结果示意图。

图8为短波候选槽点筛选示意图。

图9为涡度筛选示意图。

图10为聚类结果示意图。

图11为排除标记区域的候选槽点示意图。

图12为对各生成类进行曲线拟合示意图。

图13为拓扑类型特征图。

图14为风向弧度坐标系示意图。

图15为一个风场网格数据示意图。

图16为弧度偏转筛选示意图。

图17为标记长波槽周围区域示意图。

具体实施方式

本发明公开一种基于风场网格数据的槽线自动识别方法，包括以下步骤：

步骤1，长波槽生成：对风场网格数据进行特征点提取，判别逆向环流中心，同时提取候选长波槽点，对于每一个逆向环流中心，将其作为起点进行槽点追踪，输出单条长波槽，并将已使用槽点标注；

步骤2，短波槽生成：标记步骤1中长波槽区域，对风场网格数据进行涡度、角度偏转筛选得到候选短波槽点，排除长波槽区域内的候选短波槽点后，对剩下的候选短波槽点的集合进行最小生成树聚类，在各个生成类的候选短波槽点集合中进行曲线拟合作为短波槽分析结果；

步骤3，计算机气象软件根据步骤2的结果，自动绘制风场的槽线并显示。

其中，步骤1包括以下步骤：

步骤1-1，提取逆向环流中心：在风场网格数据中查找局部风速最小点得到特征点，计算特征点处雅可比Jacobian矩阵并判断矢量拓扑结构，得到逆向环流中心；

步骤1-2，提取候选长波槽点：遍历风场网格数据，计算网格点风矢方向，依据水平、垂直两个方向上相邻网格风矢方向特点判断候选槽点，得到候选槽点集合；

步骤1-3，追踪槽线：通过步骤1-1、步骤1-2分别得到逆向环流中心和候选槽点集合，将用于南向追踪的候选槽点记为WP，WP构成的候选槽点集合为S_W，将用于东向追踪的候选槽点记为EP，EP构成的候选槽点集合为S_E，由长波槽起点开始分两个方向分别追踪两类候选槽点，分别得到槽线分析结果为点列PList_W和PList_E，即单条长波槽。

步骤1-1包括以下步骤：

步骤1-1-1，查找特征点：遍历风场网格数据，在每一个3*3网格中心查找局部风速最小点，选取对称位置上风向相反的点作为特征点；

步骤1-1-2，筛选逆向环流中心：构建Jacobian矩阵，求矩阵特征值，通过特征值判断特征点类型，筛选出所需逆向环流中心。

步骤1-1-1包括以下步骤：

步骤1-1-1-1，查找局部风速最小点：从风场网格数据左上角开始选取一个网格数据3*3点阵，寻找局部风速最小点，该点阵形式如下表所示：

V11	V12	V13
			V21	V22	V23
V31	V32	V33

每个方格代表风场中的一个网格点，V₁₁,V₁₂,...V₃₃代表在各自网格点的风速大小，由下式计算：

$V_{ii,jj}=\sqrt{u_{ii,jj}^2+v_{ii,jj}^2},$

其中，ii、jj分别为网格点行、列坐标索引，ii、jj取值范围1～3，u_ij,v_ij分别为网格点坐标(ii,jj)处的风速水平分量与垂直分量。在该点阵中，如果中心点V₂₂的值比它周围的8个点处的值都小，则判别该点为局部风速最小点，并将其作为一个特征点候选点；

步骤1-1-1-2，结合风向判断：对于满足步骤1-1-1-1的特征点候选点，判断关于该特征点候选点成中心对称的网格点上的风速方向是否相反，即在该3*3矢量点阵中，是否满足以下条件：

u₁₁×u₃₃＜0且v₁₁×v₃₃＜0，

或

u₁₃×u₃₁＜0且v₁₃×v₃₁＜0，

如满足，该候选点即为特征点，否则判定不是特征点；

步骤1-1-1-3，按行列顺序坐标增1对风场网格数据所有3*3点阵，重复步骤1-1-1-1～1-1-1-2直到最后一个网格结束；

步骤1-1-2包括以下步骤：

步骤1-1-2-1，求解Jacobian矩阵J：

$J=\left(\begin{array}{cc}\frac{\∂u}{\∂x}& \frac{\∂u}{\∂y}\\ \frac{\∂v}{\∂x}& \frac{\∂v}{\∂y}\end{array}\right),$

其中u,v分别为某一点处风速水平分量与垂直分量，分别表示该点风速水平分量u在水平与垂直方向上的偏导；分别表示该点风速垂直分量v在水平与垂直方向上的偏导。

对于矩阵中非边界数据点，即下标大于1，行标小于矩阵行数，列表小于矩阵列数的网格点，坐标索引为(i,j)的网格点处偏导用下式表示：

$\frac{\∂u_{i,j}}{\∂x}=\frac{1}{2}(u_{i,j+1}-u_{i,j-1})\frac{\∂u_{i,j}}{\∂y}=\frac{1}{2}(u_{i-1,j}-u_{i+1,j}),$

$\frac{\∂v_{i,j}}{\∂x}=\frac{1}{2}(v_{i,j+1}-v_{i,j-1})\frac{\∂v_{i,j}}{\∂y}=\frac{1}{2}(v_{i-1,j}-v_{i+1,j}),$

首行首列点处的偏导求解采取下式表示：

$\frac{\∂u_{i,1}}{\∂x}=u_{i,2}-u_{i,1}\frac{\∂u_{1,j}}{\∂y}=u_{1,j}-u_{2,j},$

$\frac{\∂v_{i,1}}{\∂x}=v_{i,2}-v_{i,1}\frac{\∂v_{1,j}}{\∂y}=v_{1,j}-v_{2,j},$

末行Row、末列Col点处的偏导求解采用下式表示：

$\frac{\∂u_{i,Col}}{\∂x}=u_{i,Col}-u_{i,Col-1}\frac{\∂u_{Row,j}}{\∂y}=u_{Row-1,j}-u_{Row,j},$

$\frac{\∂v_{i,Row}}{\∂x}=v_{i,Row}-v_{i,Row-1}\frac{\∂v_{Col,j}}{\∂y}=v_{Col-1,j}-v_{Col,j},$

步骤1-1-2-2，求解特征值：采用下式求解坐标索引为(i,j)的网格点处Jacobian矩阵J_i,j：

$J_{i,j}=\left(\begin{array}{cc}\frac{\∂u_{i,j}}{\∂x}& \frac{\∂u_{i,j}}{\∂y}\\ \frac{\∂v_{i,j}}{\∂x}& \frac{\∂v_{i,j}}{\∂y}\end{array}\right),$

求解行列式：

$\left|\begin{array}{cc}\frac{\∂u_{i,j}}{\∂x}-\mathrm{\λ}& \frac{\∂u_{i,j}}{\∂y}\\ \frac{\∂v_{i,j}}{\∂x}& \frac{\∂v_{i,j}}{\∂y}-\mathrm{\λ}\end{array}\right|=0,$

得到该点的Jacobian矩阵的特征值，即λ的两个解，记为λ₁,λ₂；

步骤1-1-2-3，判断特征点类型：根据矢量拓扑结构特性，提取其中吸引旋转点和旋转点两种拓扑结构即为逆向环流中心，其特征形式如图13所示：

通过步骤1-1-2-1计算得到的特征值λ₁,λ₂，其中λ₁对应的实部与虚部分别为R₁,I₁，λ₂对应实部和虚部分别为R₂,I₂，当R₁＜0,R₂＜0；I₁≠0,I₂≠0时，判定特征点为吸引旋转点；当R₁＝0,R₂＝0；I₁≠0,I₂≠0时，判定特征点为旋转点；

步骤1-1-2-4，对于每一个特征点，重复步骤1-1-2-1～1-1-2-3直到最后一个特征点，得到逆向环流中心点集合。

步骤1-2包括以下步骤：

步骤1-2-1，计算风矢方向：在风场网格数据中，建立如图14所示的弧度坐标系，北风弧度为9东°，东风的弧度为180°，南风弧度为270°。通过风矢水平、垂直两个方向上的分量求得每一个网格点处风矢的方向。

风场网格数据的网格尺寸为Row*Col，P_i,j表示坐标索引为(i,j)的网格点，P_i,j.X为P_i,j点处X方向上的几何坐标，P_i,j.Y为P_i,j点处Y方向上的几何坐标，P_i,j处风矢方向为α_i,j，分别计算一个由四个相邻网格点组成的网格中各点处弧度，即在图15所示网格中，计算P_i,j、P_i,j+1、P_i+1,j、P_i+1,j+1处的方向为α_i,j、α_i,j+1、α_i+1,j、α_i+1,j+1；

步骤1-2-2，判断α_i,j与α_i,j+1角度范围：若α_i,j＜90°且270°＜α_i,j+1＜360°，则P_i,j、P_i,j+1两点间存在候选槽点WP，候选槽点WP在X方向上、Y方向上的几何坐标WP.X和WP.Y分别为：

$WP.X=\frac{180-{\mathrm{\α}}_{i,j}}{{\mathrm{\α}}_{i,j+1}-{\mathrm{\α}}_{i,j}}*(P_{i,j+1}.X-P_{i,j}.X)+P_{i,j}.X,$

WP.Y＝P_i,j.Y；

步骤1-2-3，判断α_i,j与α_i+1,j角度范围：若180°＜α_i,j＜270°且270°＜α_i+1,j＜360°，那么P_i,j、P_i+1,j两点间存在候选槽点EP，候选槽点EP在X方向上、Y方向上的几何坐标EP.X和EP.Y分别为：

EP.X＝P_i,j.X，

$EP.Y=\frac{180-{\mathrm{\α}}_{i,j}}{{\mathrm{\α}}_{i+1,j}-{\mathrm{\α}}_{i,j}}*(P_{i+1,j}.Y-P_{i,j}.Y)+P_{i,j}.Y;$

步骤1-2-4，遍历风场网格中所有网格点数据，重复步骤1-2-1～1-2-3筛选候选长波槽点直到最后一个网格，得到候选槽点集合。

步骤1-3包括以下步骤：

步骤1-3-1，选择追踪起点：从逆向环流中心集合中按序任意选取一个点作为长波槽追踪起点StartP，当追踪方向Direction＝270°时，设候选槽点集合S＝S_W，点列PList＝PList_W；当追踪方向Direction＝180°，设候选槽点集合S＝S_E，点列PList＝PList_E；

步骤1-3-2，筛选候选长波槽点：通过距离、方向筛选选取符合条件的候选槽点：

步骤1-3-2-1，距离筛选：遍历候选槽点集合S，通过距离阈值ThredL对其进行筛选，设水平相邻网格点间距为disAdjacent，网格点间距是经纬度上网格点风矢对应到绘制几何坐标系中两个水平相邻网格点间的几何距离，由绘制风矢时选取的距离决定。ThredL应为3*disAdjacent～5*disAdjacent之间。点P_k为S集合中第k点，点P_k与追踪起点StartP距离为L_k，若L_k＜ThredL，将P_k加入筛选后候选槽点集合S1中，否则舍弃，若不存在满足L_k＜ThredL的点，返回执行步骤1-3-1；

步骤1-3-2-2，方向筛选：遍历候选槽点集合S1，点P_m为候选槽点集合S1中第m点，点P_m与追踪起点StartP构成的直线与追踪方向Direction夹角为γ_m，

为了保证槽线是向正确的方向延伸，同时保证槽线平滑特性，候选槽点与追踪方向构成的弧度必须满足弧度偏转的要求。如图16所示，若槽线已探索延伸至点P₁，追踪方向为AB，下一个可能在槽线上的关键点应在以AB为轴一定大小的张角范围内。那么左侧存在一个可偏转的弧度，其阈值设定为α_left，右侧存在一个可偏转弧度阈值，其阈值设定为α_right，设置α_left＝30°，α_right＝30°。槽线的追踪设置在上述规定角度范围内。

将满足上述条件且γ_m＜Thredγ的P_m加入筛选后的候选槽点集合S2，若不存在满足上述条件的点，返回步骤1-3-1；

步骤1-3-3，通过评估方程确定下一候选点：设置评估方程Evaluation(L,α)对每一个符合条件的候选槽点进行评价：

$Evaluation(L,{\mathrm{\α}}_{field})=\frac{WeightsA}{L^2}+\frac{WeightsB}{|{\mathrm{\α}}_{field}-{\mathrm{\α}}_{fixed}|},$

其中，L为候选槽点与已延伸槽线上最后一点的距离，α_field为候选槽点与最后一点构成线段与追踪方向的夹角，α_fixed为使得槽线追踪向着需要的方向偏转的固定偏转量，WeightsA和WeightsB分别为距离判断和弧度偏转判断权值，这两个权值的设定是对距离和弧度偏转两个要素的重要性进行权衡，取值范围均为0～1。遍历候选槽点集合S2，查找使得Evaluation(L,α_field)最大的点P_Next，将点P_Next加入点列PList；

步骤1-3-4，重确定追踪起点、追踪方向：求取追踪起点StartP与下一个点P_Next在坐标系中角度δ，设置追踪方向Direction＝δ+α_fixed，将追踪起点StartP设置为P_Next；

步骤1-3-5，按序遍历逆向环流中心集合，重复步骤1-3-2、1-3-3、1-3-4直到最后一个逆向环流中心点。

步骤2包括以下步骤：

步骤2-1，标记长波槽区域：由步骤1得到组成长波槽的点列，对这些点周围一定范围进行标记，该范围距离应该保证在4*disAdjacent～6*disAdjacent之间，disAdjacent由步骤1-3-2-1得到。得到标记的区域M；如图17所示，OA、OB是长波槽分析结果，标记A₁O₁B₁与A₂O₂B₂之间的区域；

步骤2-2，筛选候选短波槽点：从风场网格数据中按角度偏转和涡度筛选阈值查找大于阈值的点，即风切变较大的候选短波槽点；

步骤2-2-1，筛选偏转角度：对风场网格数据中的每一个网格点，分别在水平和垂直两个方向上求取相邻网格点风矢的方向偏转大小；

步骤2-2-1-1，求取角度偏转量：Angle_horizon_s,t、Angle_vertical_s,t分别为网格点(s,t)处水平与垂直方向上的角度偏转，则：

Angle_horizon_s,t＝α_s,t-α_s,t+1，

Angle_vertical_s,t＝α_s,t-α_s+1,t，

若Angle_horizon_s,t＜0，设置Angle_horizon_s,t＝-Angle_horizon_s,t；

若Angle_vertical_s,t＜0，设置Angle_vertical_s,t＝-Angle_vertical_s,t；

若Angle_vertical_s,t＞180°，设置Angle_vertical_s,t＝360°-Angle_vertical_s,t；

若Angle_horizon_s,t＞180°，设置Angle_horizon_s,t＝360°-Angle_horizon_s,t；

步骤2-2-1-2，完成步骤2-2-1-1角度偏转计算后，计算偏转量的筛选阈值AngleThred：

$\begin{array}{c}AngleThred=\mathrm{\β}(\frac{1}{Col*(Row-1)}\underset{s,t}{\Σ}Angle\_{\mathrm{horizon}}_{s,t}\\ +\frac{1}{(Col-1)*Row}\underset{s,t}{\Σ}Angle\_{\mathrm{vertical}}_{s,t})\end{array},$

其中β为可调整权值，取值范围为0～5；

步骤2-2-1-3，遍历所有网格点，若Angle_vertical_s,t＞AngleThred，则将网格中第s行、第t列点与第s行、第t+1列点之间中点标记为候选槽点，加入候选短波槽点集合O₁；若Angle_horizon_s,t＞AngleThred，则将网格中第s行、第t列点与第s+1行、第t列点之间中点标记为候选槽点，加入候选短波槽点集合O₁；

步骤2-2-2，筛选涡度：遍历风场网格中网格点数据求得各点处涡度值，计算涡度筛选阈值进行筛选提取候选短波槽点：

步骤2-2-2-1，在风场网格数据中，点(i,j)处涡度大小Vorticity_ij为：

Vorticity_i,j＝|u_x-v_y|

其中u_x是水平分量u在水平方向上的偏导，v_y是垂直分量v在垂直方向上的偏导，同步骤1-1-2-1用差分代替上式中的偏微分计算涡度大小；

步骤2-2-2-2，计算涡度筛选阈值VorticityThred：

$VorticityThred=\mathrm{\γ}*\left(\frac{1}{Row*Col}\underset{i,j}{\Σ}{\mathrm{Vorticity}}_{i,j}\right),$

其中，γ是可调整的涡度筛选权重，取值范围为0～5；

步骤2-2-2-3，标记候选短波槽点：若Vorticity_i,j＞VorticityThred，则将风场网格中第i行、第j列点标记为候选短波槽点，将其加入候选短波槽点集合O₁；

步骤2-3，排除部分候选短波槽点：结合步骤2-1已标记的部分区域M，筛选后的候选槽点集合为O₂：

O₂＝O₁-(O₁∩M)；

步骤2-4，候选槽点最小生成树聚类：对候选槽点集合O₂中的点依据点集合坐标位置建立加权无向连通图，进行最小生成树聚类，生成多个分类：

步骤2-5，曲线拟合：对生成的分类设定拟合多项式，通过最小二乘法确定参数，实现多项式曲线拟合，提取短波槽分析结果。

步骤2-4包括以下步骤：

步骤2-4-1，建立加权无向连通图G＝(Q,E)：将候选槽点集合O₂中的槽点按需排列编号为q₁,q₂,...,q_n，n为候选槽点集合O₂中点的个数，顶点集合为Q＝{q₁,q₂,...,q_n}，边集合为E，E＝{＜q₁,q₂＞,＜q₁,q₃＞,...,＜q_n-1,q_n＞}，其中＜q_i,q_j＞为顶点q_i与q_j的距离；

步骤2-4-2，初始化顶点集合：初始化顶点集合Q_new＝{q_z}，其中q_z为Q中任意一点，初始边集合E_new＝{}为空；

步骤2-4-3，构建最小生成树：从初始化顶点集合Q_new＝{q_z}中的起点开始，在边集合E中查找与顶点集合Q_new中的点构成的边中权值最小的边，如果这条边上的点不在顶点集合Q_new中，则将该点加入到顶点集合Q_new中：

步骤2-4-3-1，在边集合E中选取权值最小的边＜q_l,q_d＞，其中q_l为顶点集合Q_new中的元素，q_d不在顶点集合Q_new中，如果存在有多条边权值均为最小值，则任意选取其中之一；

步骤2-4-3-2，将q_d加入到顶点集合Q_new中，将＜q_l,q_d＞边加入到边集合E_new中；

步骤2-4-3-3，重复步骤2-4-3-1～2-4-3-2直到Q_new＝Q，输出聚类结果，使用集合Q_new、E_new描述所得最小生成树；

步骤2-4-4，分割子树：对步骤2-4-3-3得到的最小生成树边集合E_new，设置筛选阈值ClustThred，，ClustThred应在3*disAdjacent～5*disAdjacent范围内。遍历边集合E_new所有元素，若其权值大于ClustThred，则去掉该条边，将最小生成树分为两个以上的分树，即可得到候选短波槽点集合通过最小生成树方法划分的子类。

步骤2-5包括以下步骤：

步骤2-5-1，设定多项式：依据多项式f(h)＝a₁h⁴+a₂h³+a₃h²+a₄h+a₅进行拟合，其中h为自变量，f(h)为因变量，a₁,a₂,a₃,a₄,a₅为待确定参数，对一个总共g个点的分类点集合Sub＝{p₁,p₂,...,p_g}，每一个点由水平方向上和垂直方向上两个坐标决定，分别设为p_r.X,p_r.Y(r＝1,2,...,g)，将这些点水平方向上的坐标作为自变量，垂直方向上的坐标作为因变量；

步骤2-5-2，采用最小二乘法求取参数：依据偏差平方和最小的原则求取待确定参数，

距离平方和即

分别对等式右侧a₁,a₂,a₃,a₄,a₅求偏导得：

令上式都为零得到范德蒙德矩阵：

$\left[\begin{array}{ccccc}g& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^7\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^7& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^8\end{array}\right]*\left[\begin{array}{c}{a}_5\\ a_4\\ a_3\\ a_2\\ a_1\end{array}\right]=\left[\begin{array}{c}\underset{r=1}{\overset{g}{\Σ}}{p}_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4*p_r.Y\end{array}\right],$

将该范德蒙德矩阵简化后得到：

$\left[\begin{array}{ccccc}1& p_1.X& p_1.X^2& p_1.X^3& p_1.X^4\\ 1& p_2.X& p_2.X^2& p_2.X^3& p_2.X^4\\ 1& p_3.X& p_3.X^2& p_3.X^3& p_3.X^4\\ .& .& .& .& .\\ .& .& .& .& .\\ 1& p_g.X& p_g.X^2& p_g.X^3& p_g.X^4\end{array}\right]*\left[\begin{array}{c}{a}_5\\ a_4\\ a_3\\ a_2\\ a_1\end{array}\right]=\left[\begin{array}{c}{p}_1.Y\\ p_2.Y\\ p_3.Y\\ .\\ .\\ p_g.Y\end{array}\right],$

求得a₁,a₂,a₃,a₄,a₅的值，即得多项式拟合曲线，将其作为各分类点集合短波槽分析结果。

实施例1

本实施例包括长波槽提取和短波槽提取，其处理流程图如图1所示，整个方法分为长波槽、短波槽提取两个主要步骤，下面分别介绍实施例部分的主要流程。

1.长波槽提取

首先对风场数据进行提取，截取N 0°---N 60°，E60°---E180°之间包含我国区域的数据绘制出风场图，如图4a、图4b所示。图4a为2006年5月1日00时500hpa风矢绘制结果。图4b中为2006年5月2日00时500hpa风矢绘制结果。选取图4a中2006年5月1日00时500hpa数据进行进一步的处理。

按照步骤11提取逆向环流中心，如图5所示，黑色点为提取得到的逆向环流中心。

按照步骤12提取长波候选槽点，如图6所示，分别判断得出东向追踪候选槽点、南向追踪候选槽点。

按照步骤13对长波槽进行追踪得到追踪结果，如图7所示。

2.短波槽提取

长波槽提取结束后，首先对长波槽区域进行标记，再按照步骤22进行短波候选槽点筛选，如图8所示，分别求出水平方向上风矢方向变化大的点、垂直方向上风矢变化大的点，并标识出来。图9为对涡度的筛选，由于数据截取北半球部分，因此逆时针旋转涡度为正，顺时针旋转涡度为负，通过阈值筛选出涡度大的点。

按照步骤24，对候选槽点进行最小生成树聚类，其聚类结果如图10所示。通过步骤21的筛选，去掉标记区域的候选槽点，如图11所示，R₁,R₂,R₃,R₄为需要排除区域。

按照步骤25，对各生成类进行曲线拟合，如图12所示。

Claims (10)

1.一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，包括以下步骤：

步骤1，长波槽生成：对风场网格数据进行特征点提取，判别逆向环流中心，同时提取候选长波槽点，对于每一个逆向环流中心，将其作为起点进行槽点追踪，输出单条长波槽，并将已使用槽点标注；

步骤2，短波槽生成：标记步骤1中长波槽区域，对风场网格数据进行涡度、角度偏转筛选得到候选短波槽点，排除长波槽区域内的候选短波槽点后，对剩下的候选短波槽点的集合进行最小生成树聚类，在各个生成类的候选短波槽点集合中进行曲线拟合作为短波槽分析结果；

步骤3，计算机气象软件根据步骤2的结果，自动绘制风场的槽线并显示。

2.根据权利要求1所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤1包括以下步骤：

步骤1-1，提取逆向环流中心：在风场网格数据中查找局部风速最小点得到特征点，计算特征点处雅可比Jacobian矩阵并判断矢量拓扑结构，得到逆向环流中心；

步骤1-2，提取候选长波槽点：遍历风场网格数据，计算网格点风矢方向，依据水平、垂直两个方向上相邻网格风矢方向判断候选槽点，得到候选槽点集合；

步骤1-3，追踪槽线：将用于南向追踪的候选槽点记为WP，WP构成的候选槽点集合为S_W，将用于东向追踪的候选槽点记为EP，EP构成的候选槽点集合为S_E，由长波槽起点开始分两个方向分别追踪两类候选槽点，分别得到槽线分析结果为点列PList_W和点列PList_E，即单条长波槽。

3.根据权利要求2所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤1-1包括以下步骤：

步骤1-1-1，查找特征点：遍历风场网格数据，在每一个3*3网格中心查找局部风速最小点，选取对称位置上风向相反的点作为特征点；

步骤1-1-2，筛选逆向环流中心：构建Jacobian矩阵，求矩阵特征值，通过特征值判断特征点类型，筛选出所需逆向环流中心。

4.根据权利要求3所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤1-1-1包括以下步骤：

步骤1-1-1-1，查找局部风速最小点：从风场网格数据左上角开始选取一个网格数据3*3点阵，寻找局部风速最小点，该点阵形式如下表所示：

V₁₁ V₁₂ V₁₃ V₂₁ V₂₂ V₂₃ V₃₁ V₃₂ V₃₃

每个方格代表风场中的一个网格点，V₁₁,V₁₂,...V₃₃代表在各自网格点的风速大小，由下式计算：

$V_{ii,jj}=\sqrt{u_{ii,jj}^2+v_{ii,jj}^2},$

其中，ii、jj分别为网格点行、列坐标索引，ii、jj取值范围1～3，u_ii,jj,v_ii,jj分别为网格点坐标(ii,jj)处的风速水平分量与垂直分量，在该点阵中，如果中心点V₂₂的值比它周围的8个点处的值都小，则判别该点为局部风速最小点，并将其作为一个特征点候选点；

步骤1-1-1-2，结合风向判断：对于满足步骤1-1-1-1的特征点候选点，判断关于该特征点候选点成中心对称的网格点上的风速方向是否相反，即在该3*3矢量点阵中，是否满足以下条件：

u₁₁×u₃₃＜0且v₁₁×v₃₃＜0，

或

u₁₃×u₃₁＜0且v₁₃×v₃₁＜0，

如满足，该候选点即为特征点，否则判定不是特征点；

步骤1-1-1-3，按行列顺序坐标增1对风场网格数据所有3*3点阵，重复步骤1-1-1-1～1-1-1-2直到最后一个网格结束。

5.根据权利要求4所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤1-1-2包括以下步骤：

步骤1-1-2-1，求解Jacobian矩阵J：

$J=\left(\begin{array}{cc}\frac{\∂u}{\∂x}& \frac{\∂u}{\∂y}\\ \frac{\∂v}{\∂x}& \frac{\∂v}{\∂y}\end{array}\right),$

其中u,v分别为任意一点处风速水平分量与垂直分量，分别表示该点风速水平分量u在水平与垂直方向上的偏导；分别表示该点风速垂直分量v在水平与垂直方向上的偏导，

对于矩阵中非边界数据点，即下标大于1，行标小于矩阵行数，列标小于矩阵列数的网格点，坐标索引为(i,j)的网格点处偏导用下式表示：

$\frac{\∂u_{i,j}}{\∂x}=\frac{1}{2}(u_{i,j+1}-u_{i,j-1}),\frac{\∂u_{i,j}}{\∂y}=\frac{1}{2}(u_{i-1,j}-u_{i+1,j}),$

$\frac{\∂v_{i,j}}{\∂x}=\frac{1}{2}(v_{i,j+1}-v_{i,j-1}),\frac{\∂v_{i,j}}{\∂y}=\frac{1}{2}(v_{i-1,j}-v_{i+1,j}),$

矩阵首行首列点处的偏导求解采取下式表示：

$\frac{\∂u_{i,1}}{\∂x}=u_{i,2}-u_{i,1},\frac{\∂u_{1,j}}{\∂y}=u_{1,j}-u_{2,j},$

$\frac{\∂v_{i,1}}{\∂x}=v_{i,2}-v_{i,1},\frac{\∂v_{1,j}}{\∂y}=v_{1,j}-v_{2,j},$

矩阵末行Row、末列Col点处的偏导求解采用下式表示：

$\frac{\∂u_{i,Col}}{\∂x}=u_{i,Col}-u_{i,Col-1},\frac{\∂u_{Row,j}}{\∂y}=u_{Row-1,j}-u_{Row,j},$

$\frac{\∂v_{i,Row}}{\∂x}=v_{i,Row}-v_{i,Row-1},\frac{\∂v_{Col,j}}{\∂y}=v_{Co1-1,j}-v_{Col,j};$

步骤1-1-2-2，求解特征值：采用下式求解坐标索引为(i,j)的网格点处Jacobian矩阵J_i,j：

$J_{i,j}=\left(\begin{array}{cc}\frac{\∂u_{i,j}}{\∂x}& \frac{\∂u_{i,j}}{\∂y}\\ \frac{\∂v_{i,j}}{\∂x}& \frac{\∂v_{i,j}}{\∂y}\end{array}\right),$

求解行列式：

$\left|\begin{array}{cc}\frac{\∂u_{i,j}}{\∂x}-\mathrm{\λ}& \frac{\∂u_{i,j}}{\∂y}\\ \frac{\∂v_{i,j}}{\∂x}& \frac{\∂v_{i,j}}{\∂y}-\mathrm{\λ}\end{array}\right|=0,$

得到该点的Jacobian矩阵的特征值，即λ的两个解，记为λ₁,λ₂；

步骤1-1-2-3，判断特征点类型：根据矢量拓扑结构特性，提取其中吸引旋转点和旋转点两种拓扑结构即为逆向环流中心；

通过步骤1-1-2-1计算得到的特征值λ₁,λ₂，其中λ₁对应的实部与虚部分别为R₁,I₁，λ₂对应实部和虚部分别为R₂,I₂，当R₁＜0,R₂＜0；I₁≠0,I₂≠0时，判定特征点为吸引旋转点；当R₁＝0,R₂＝0；I₁≠0,I₂≠0时，判定特征点为旋转点；

步骤1-1-2-4，对于每一个特征点，重复步骤1-1-2-1～1-1-2-3直到最后一个特征点，得到逆向环流中心点集合。

6.根据权利要求5所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤1-2包括以下步骤：

步骤1-2-1，计算风矢方向：在风场网格数据中，建立弧度坐标系，北风弧度为90°，东风的弧度为180°，南风弧度为270°，通过风矢水平、垂直两个方向上的分量求得每一个网格点处风矢的方向，风场网格数据的网格尺寸为Row*Col，P_i,j表示坐标索引为(i,j)的网格点，P_i,j.X为P_i,j点处X方向上的几何坐标，P_i,j.Y为P_i,j点处Y方向上的几何坐标，P_i,j处风矢方向为α_i,j，分别计算一个由四个相邻网格点组成的网格中各点处弧度，即计算P_i,j、P_i,j+1、P_i+1,j、P_i+1,j+1处的方向α_i,j、α_i,j+1、α_i+1,j、α_i+1,j+1；

步骤1-2-2，判断α_i,j与α_i,j+1角度范围：若α_i,j＜90°且270°＜α_i,j+1＜360°，则P_i,j、P_i,j+1两点间存在用于南向追踪的候选槽点WP，候选槽点WP在X方向上、Y方向上的几何坐标WP.X和WP.Y分别为：

$WP.X=\frac{180-{\mathrm{\α}}_{i,j}}{{\mathrm{\α}}_{i,j+1}-{\mathrm{\α}}_{i,j}}*(P_{i,j+1}.X-P_{i,j}.X)+P_{i,j}.X,$

WP.Y＝P_i,j.Y；

步骤1-2-3，判断α_i,j与α_i+1,j角度范围：若180°＜α_i,j＜270°且270°＜α_i+1,j＜360°，那么P_i,j、P_i+1,j两点间存在用于东向追踪的候选槽点EP，候选槽点EP在X方向上、Y方向上的几何坐标EP.X和EP.Y分别为：

EP.X＝P_i,j.X，

$EP.Y=\frac{180-{\mathrm{\α}}_{i,j}}{{\mathrm{\α}}_{i+1,j}-{\mathrm{\α}}_{i,j}}*(P_{i+1,j}.Y-P_{i,j}.Y)+P_{i,j}.Y;$

步骤1-2-4，遍历风场网格中所有网格点数据，重复步骤1-2-1～1-2-3筛选候选长波槽点直到最后一个网格，得到候选槽点集合。

7.根据权利要求6所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤1-3包括以下步骤：

步骤1-3-1，选择追踪起点：从逆向环流中心集合中按序任意选取一个点作为长波槽追踪起点StartP，当追踪方向Direction＝270°时，设候选槽点集合S＝S_W，点列PList＝PList_W；当追踪方向Direction＝180°，设候选槽点集合S＝S_E，点列PList＝PList_E；

步骤1-3-2，筛选候选长波槽点：通过距离、方向筛选选取符合条件的候选槽点：

步骤1-3-2-1，距离筛选：遍历候选槽点集合S，通过距离阈值ThredL对其进行筛选，设水平相邻网格点间距为disAdjacent，阈值ThredL的值在3*disAdjacent～5*disAdjacent之间，点P_k为S集合中第k点，点P_k与追踪起点StartP距离为L_k，若L_k＜ThredL，将P_k加入筛选后候选槽点集合S1中，否则舍弃，若不存在满足L_k＜ThredL的点，返回执行步骤1-3-1；

步骤1-3-2-2，方向筛选：遍历候选槽点集合S1，点P_m为候选槽点集合S1中第m点，点P_m与追踪起点StartP构成的直线与追踪方向Direction夹角为γ_m，

候选槽点与追踪方向构成的弧度满足以下弧度偏转的条件：

若槽线已探索延伸至点P₁，追踪方向为AB，下一个可能在槽线上的关键点在以追踪方向AB为轴一定大小的张角范围内，那么左侧存在一个可偏转的弧度，其阈值设定为α_left，右侧存在一个可偏转弧度，其阈值设定为α_right，设置α_left＝30°，α_right＝30°，槽线的追踪设置在上述规定角度范围内；

将满足上述条件且γ_m＜Thredγ的点P_m加入筛选后的候选槽点集合S2，若不存在满足上述条件的点，返回步骤1-3-1；

步骤1-3-3，通过评估方程确定下一候选点：设置评估方程Evaluation(L,α)对每一个符合条件的候选槽点进行评价：

$Evaluation\left(L,{\mathrm{\α}}_{field}\right)=\frac{WeightsA}{L^2}+\frac{WeightsB}{\left|{\mathrm{\α}}_{field}-{\mathrm{\α}}_{fixed}\right|},$

其中，L为候选槽点与已延伸槽线上最后一点的距离，α_field为候选槽点与最后一点构成线段与追踪方向的夹角，α_fixed为使得槽线追踪向着需要的方向偏转的固定偏转量，WeightsA和WeightsB分别为距离判断和弧度偏转判断权值，取值范围均为0～1，遍历候选槽点集合S2，查找使得Evaluation(L,α_field)最大的点P_Next，将点P_Next加入点列PList；

步骤1-3-4，重确定追踪起点、追踪方向：求取追踪起点StartP与下一个点P_Next在坐标系中角度δ，设置追踪方向Direction＝δ+α_fixed，将追踪起点StartP设置为P_Next；

步骤1-3-5，按序遍历逆向环流中心集合，重复步骤1-3-2～1-3-4直到最后一个逆向环流中心点。

8.根据权利要求7所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤2包括以下步骤：

步骤2-1，标记长波槽区域：由步骤1-3得到组成长波槽的点列，对这些点周围范围进行标记，该范围距离限定在4*disAdjacent～6*disAdjacent之间，disAdjacent由步骤1-3-2-1得到，得到标记的区域M；

步骤2-2，筛选候选短波槽点，包括如下步骤：

步骤2-2-1，筛选偏转角度：对风场网格数据中的每一个网格点，分别在水平和垂直两个方向上求取相邻网格点风矢的方向偏转大小，包括如下步骤：

步骤2-2-1-1，求取角度偏转量：Angle_horizon_s,t、Angle_vertical_s,t分别为网格点(s,t)处水平与垂直方向上的角度偏转，则：

Angle_horizon_s,t＝α_s,t-α_s,t+1，

Angle_vertical_s,t＝α_s,t-α_s+1,t，

若Angle_horizon_s,t＜0，设置Angle_horizon_s,t＝-Angle_horizon_s,t；

若Angle_vertical_s,t＜0，设置Angle_vertical_s,t＝-Angle_vertical_s,t；

若Angle_vertical_s,t＞180°，设置Angle_vertical_s,t＝360°-Angle_vertical_s,t；

若Angle_horizon_s,t＞180°，设置Angle_horizon_s,t＝360°-Angle_horizon_s,t；

步骤2-2-1-2，完成步骤2-2-1-1角度偏转计算后，计算偏转量的筛选阈值AngleThred：

$\begin{array}{c}AngleThred=\mathrm{\β}(\frac{1}{Col*\left(Row-1\right)}\underset{s,t}{\Σ}Angle\_{\mathrm{horizon}}_{s,t}\\ +\frac{1}{\left(Col-1\right)*Row}\underset{s,t}{\Σ}Angle\_{\mathrm{vertical}}_{s,t})\end{array},$

其中β为可调整权值，取值范围为0～5；

步骤2-2-1-3，遍历所有网格点，若Angle_vertical_s,t＞AngleThred，则将网格中第s行、第t列点与第s行、第t+1列点之间中点标记为候选槽点，加入候选短波槽点集合O₁；若Angle_horizon_s,t＞AngleThred，则将网格中第s行、第t列点与第s+1行、第t列点之间中点标记为候选槽点，加入候选短波槽点集合O₁；

步骤2-2-2，筛选涡度：遍历风场网格中网格点数据求得各点处涡度值，计算涡度筛选阈值进行筛选提取候选短波槽点：

步骤2-2-2-1，在风场网格数据中，网格点(i,j)处涡度大小Vorticity_i,j为：

Vorticity_i,j＝|u_x-v_y|，

其中u_x是水平分量u在水平方向上的偏导，v_y是垂直分量v在垂直方向上的偏导，用差分代替上式中的偏微分计算涡度大小；

步骤2-2-2-2，计算涡度筛选阈值VorticityThred：

$VorticityThred=\mathrm{\γ}*\left(\frac{1}{Row*Col}\underset{i,j}{\Σ}{\mathrm{Vorticity}}_{i,j}\right),$

其中，γ是可调整的涡度筛选权重，取值范围为0～5；

步骤2-2-2-3，标记候选短波槽点：若Vorticity_i,j＞VorticityThred，则将风场网格中第i行、第j列点标记为候选短波槽点，将其加入候选短波槽点集合O₁；

步骤2-3，排除部分候选短波槽点：结合步骤2-1已标记的部分区域M，筛选后的候选槽点集合为O₂：

O₂＝O₁-(O₁∩M)；

步骤2-4，候选槽点最小生成树聚类：对候选槽点集合O₂中的点依据点集合坐标位置建立加权无向连通图，进行最小生成树聚类，生成多个分类：

步骤2-5，曲线拟合：对生成的分类设定拟合多项式，通过最小二乘法确定参数，实现多项式曲线拟合，提取短波槽分析结果。

9.根据权利要求8所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤2-4包括以下步骤：

步骤2-4-1，建立加权无向连通图G＝(Q,E)：将候选槽点集合O₂中的槽点按需排列编号为q₁,q₂,...,q_n，n为候选槽点集合O₂中点的个数，顶点集合为Q＝{q₁,q₂,...,q_n}，边集合为E，E＝{＜q₁,q₂＞,＜q₁,q₃＞,...,＜q_n-1,q_n＞}，其中＜q_i,q_j＞为顶点q_i与q_j的距离；

步骤2-4-2，初始化顶点集合：初始化顶点集合Q_new＝{q_z}，其中q_z为_Q中任意一点，z取值范围为自然数，初始边集合E_new＝{}为空；

步骤2-4-3，构建最小生成树：从初始化顶点集合Q_new＝{q_z}中的起点开始，在边集合E中查找与顶点集合Q_new中的点构成的边中权值最小的边，如果这条边上的点不在顶点集合Q_new中，则将该点加入到顶点集合Q_new中，包括：

步骤2-4-3-1，在边集合E中选取权值最小的边＜q_l,q_d＞，其中q_l为顶点集合Q_new中的元素，q_d不在顶点集合Q_new中，如果存在有多条边权值均为最小值，则任意选取其中之一；

步骤2-4-3-2，将q_d加入到顶点集合Q_new中，将＜q_l,q_d＞边加入到边集合E_new中；

步骤2-4-3-3，重复步骤2-4-3-1～2-4-3-2直到Q_new＝Q，输出聚类结果，使用集合Q_new、E_new描述所得最小生成树；

步骤2-4-4，分割子树：对步骤2-4-3-3得到的最小生成树边集合E_new，设置筛选阈值ClustThred，ClustThred限定在3*disAdjacent～5*disAdjacent范围内，遍历边集合E_new所有元素，若其权值大于ClustThred，则去掉该条边，将最小生成树分为两个以上的分树，即可得到候选短波槽点集合通过最小生成树方法划分的子类。

10.根据权利要求9所述的一种计算机气象软件中基于风场数据的槽线自动绘制方法，其特征在于，步骤2-5包括以下步骤：

步骤2-5-1，设定多项式：依据多项式f(h)＝a₁h⁴+a₂h³+a₃h²+a₄h+a₅进行拟合，其中h为自变量，f(h)为因变量，a₁,a₂,a₃,a₄,a₅为待确定参数，对一个总共g个点的分类点集合Sub＝{p₁,p₂,...,p_g}，g取值范围为自然数，每一个点由水平方向上和垂直方向上两个坐标决定，分别设为p_r.X和p_r.Y，r＝1,2,...,g，将这些点水平方向上的坐标作为自变量，垂直方向上的坐标作为因变量；

步骤2-5-2，采用最小二乘法求取参数：依据偏差平方和最小的原则求取待确定参数，

距离平方和即

分别对等式右侧a₁,a₂,a₃,a₄,a₅求偏导得：

令上式都为零得到范德蒙德矩阵：

$\left[\begin{array}{ccccc}g& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^7\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^5& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^6& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^7& \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^8\end{array}\right]*\left[\begin{array}{c}{a}_5\\ a_4\\ a_3\\ a_2\\ a_1\end{array}\right]=\left[\begin{array}{c}\underset{r=1}{\overset{g}{\Σ}}{p}_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^2*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^3*p_r.Y\\ \underset{r=1}{\overset{g}{\Σ}}{p}_r.X^4*p_r.Y\end{array}\right],$

将该范德蒙德矩阵简化后得到：

$\left[\begin{array}{ccccc}1& p_1.X& p_1.X^2& p_1.X^3& p_1.X^4\\ 1& p_2.X& p_2.X^2& p_2.X^3& p_2.X^4\\ 1& p_3.X& p_3.X^2& p_3.X^3& p_3.X^4\\ \begin{array}{c}.\\ .\end{array}& \begin{array}{c}.\\ .\end{array}& \begin{array}{c}.\\ .\end{array}& \begin{array}{c}.\\ .\end{array}& \begin{array}{c}.\\ .\end{array}\\ 1& p_g.X& p_g.X^2& p_g.X^3& p_g.X^4\end{array}\right]*\left[\begin{array}{c}{a}_5\\ a_4\\ a_3\\ a_2\\ a_1\end{array}\right]=\left[\begin{array}{c}{p}_1.Y\\ p_2.Y\\ p_3.Y\\ \begin{array}{c}.\\ .\end{array}\\ p_g.Y\end{array}\right],$

求得a₁,a₂,a₃,a₄,a₅的值，即得多项式拟合曲线，将其作为各分类点集合短波槽分析结果。