CN108875127B - 一种计算机气象软件中基于风场数据的槽线修正方法 - Google Patents

一种计算机气象软件中基于风场数据的槽线修正方法 Download PDF

Info

Publication number
CN108875127B
CN108875127B CN201810400897.5A CN201810400897A CN108875127B CN 108875127 B CN108875127 B CN 108875127B CN 201810400897 A CN201810400897 A CN 201810400897A CN 108875127 B CN108875127 B CN 108875127B
Authority
CN
China
Prior art keywords
point
points
slot
wind
slot line
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810400897.5A
Other languages
English (en)
Other versions
CN108875127A (zh
Inventor
李骞
黄岩
范茵
马烁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201810400897.5A priority Critical patent/CN108875127B/zh
Publication of CN108875127A publication Critical patent/CN108875127A/zh
Application granted granted Critical
Publication of CN108875127B publication Critical patent/CN108875127B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明公开了一种计算机气象软件中基于风场数据的槽线修正方法,包括:提取风场特征点:计算各格点的风向偏转角度及涡度,并提取风场特征点;槽线修正类型判别:分别设定槽线补充分析及槽线订正判别区域,根据判别区域中槽线和风场特征点的情况,判别槽线的修正类型;槽线补充分析:根据判别区域中风场特征点的数量,确定补充分析起点,利用最小生成树聚类算法和最小二乘法曲线拟合方法,提取补充点集并完成槽线的补充分析;槽线订正处理:确定槽线订正锚点并建立Laplacian变形框架,通过求解带位置约束的优化方程,获得变形之后各槽点的坐标,将槽点按序连接并进行平滑处理以完成槽线订正。

Description

一种计算机气象软件中基于风场数据的槽线修正方法
技术领域
本发明涉及一种计算机气象软件中基于风场数据的槽线修正方法。
背景技术
近年来,由于气象领域中预报要求的不断提高,业务部门需获取大量的观探测数据和数值模拟数据,并进行及时有效的处理,这对准确高效的分析诸多信息并完成预报工作提出极大的挑战。目前锋面、流线等特征曲线的识别绘制已能够依靠计算机自动实现,但由于高空天气系统的复杂性及诸多气象约束条件,槽线自动分析方法发展较为缓慢。
目前已有的槽线自动分析方法,大都基于气压场这一单一数据场分析提取槽线。气象槽线的形成与空气流动及等压面上的位势高低等均有关联,在实际的人工分析槽线过程中,预报员也并不会孤立的考虑单一数据场来提取槽线,而是根据风向切变和等压线走势等要素,综合分析槽线。因此,槽线自动分析的数据来源不能只局限于单一气压场,而应结合风场要素综合分析槽线,否则易导致错、漏分析槽线问题,这实际上是一个多变量数据场综合分析处理问题。
因此,针对基于单一气压场进行槽线自动分析所存在的局限性,通过融合风场要素对自动分析槽线进行修正,对于提高槽线自动分析的准确性具有重要的理论研究意义和实际应用价值。
发明内容
本发明针对现有技术不足,提供了一种计算机气象软件中基于风场数据的槽线修正方法,具体包括如下步骤:
步骤1,提取风场特征点:获取风场格点数据,根据风场格点数据计算各格点经向和纬向上的风速,由风速计算结果得到各点风向角度,并映射到角度坐标系中,根据风向角度计算各点的纬向、经向风向偏转角度以及其涡度,然后分别设定经向、纬向风向偏转角度阈值和涡度阈值,在风场格点数据中提取风场特征点;
步骤2,判别槽线修正类型:在槽线周围设定槽线补充分析判别区域,根据该区域中风场特征点的数量,判断是否需进行槽线补充分析,如需要则执行步骤3;设定槽线订正判别区域,根据该区域中风场特征点与槽线的相对位置,判断是否需进行槽线订正处理,如需要则执行步骤4;
步骤3,槽线补充分析:根据步骤2所判别得到的槽线修正类型,对于需进行补充分析的槽线,首先以补充分析判别区域对应的槽线端点作为补充分析起点,然后在风场特征点中提取槽线补充分析所需的点集,最后由补充分析起点开始,在点集中通过曲线拟合完成槽线补充分析;
步骤4,槽线订正处理:根据步骤2所判别得到的槽线修正类型,对于需进行订正的槽线,首先存储风场格点数据中所包含的网格信息(网格信息即各格点几何坐标及其之间的距离关系)并将槽线上槽点的几何坐标转换为Laplacian坐标,同时建立Laplacian坐标系,其中,Laplacian坐标的定义如下,
Figure BDA0001645681850000021
其中,δi即为顶点vi的Laplacian坐标,L()为Laplacian算子,vj为顶点vi的邻接点,N(i)为顶点vi的邻接点索引的集合,ωij为vj相对顶点vi的权值,目前常用的包括均匀权值、余切权值、正切权值等,本发明采用均匀权值,即ωij=1/di(di为vi点的度,即vi点邻接点的数量);
然后根据风场特征点和槽点确定槽线订正锚点,最后,通过求解带位置约束的优化方程,获得订正之后各槽点的几何坐标,连接槽点并平滑以完成槽线订正处理;
步骤5,根据步骤3、4的结果,输出槽线分析结果并显示在计算机屏幕上。
本发明步骤1包括以下步骤:
步骤1-1,初步提取风场特征点:于欧洲中期天气预报中心(ECMWF)官网下载获取风场格点数据,在风场格点数据中,根据其中各点的经向、纬向风速分量计算格点风向角度,并映射到角度坐标系中,然后对于全球范围内的离散风场格点数据,分别计算各格点的纬向、经向风向偏转角度,然后根据各格点的纬向、经向风向偏转角度,设置纬向、经向风向偏转角度阈值,对风场特征点进行初步提取;
步骤1-2,二次提取风场特征点:根据各格点的经向、纬向风速分量,计算格点的涡度值,然后根据格点涡度值设置涡度阈值,对风场特征点进行二次提取。
本发明步骤1-1包括以下步骤:
步骤1-1-1,通过如下公式计算格点风向角度:
Figure BDA0001645681850000022
其中,αij为第i行第j列的格点的风向角度,uij为该格点对应的纬向风速分量,vij为该格点对应的经向风速分量,由此分别计算各格点的风向角度;
步骤1-1-2,将格点风向角度映射至角度坐标系:将角度坐标系规定如下:定义正北为0°,正东为90°,正南为180°,正西为270°,根据经向、纬向的风速分量正负,将格点风向角度映射至角度坐标系中,公式如下所示:
Figure BDA0001645681850000031
映射完成之后,得到在角度坐标系下第i行第j列格点对应的风向角度值θij
步骤1-1-3,计算格点的纬向、经向风向偏转角度:在全球风场格点数据中(全球数据规模为361×720,数据分辨率为0.5×0.5,其中赤道对应的行索引为180,本初子午线对应的列索引为0),通过如下公式计算格点的纬向、经向风向偏转角度:
Figure BDA0001645681850000032
Figure BDA0001645681850000033
其中,αu_ij为第i行第j列格点的纬向风向偏转角度,αv_ij为第i行第j列格点的经向风向偏转角度,αi,1为第i行第1列格点的风向角度,αi,719为第i行第719列格点的风向角度,αi,j+1为第i行第j+1列格点的风向角度,αi+1,j为第i+1行第j列格点的风向角度;
步骤1-1-4,设置风向偏转角度阈值,对风场特征点进行初步提取:纬向、经向风向偏转角度阈值的计算公式如下所示:
Figure BDA0001645681850000034
Figure BDA0001645681850000035
其中,thresholdu、thresholdv分别为纬向、经向风向偏转角度阈值,
Figure BDA0001645681850000041
为风向偏转角度阈值系数(根据多次提取结果,在最优结果中确定系数
Figure BDA0001645681850000042
的值,一般设置为2.5),sumu和sumv分别为范围内格点纬向、经向风向偏转角度之和,row为风场格点数据行数,rank为风场格点数据列数;
规定若格点的纬向、经向风向偏转角度中,至少有一项大于其对应的风向偏转角度阈值,则将该格点初步提取为风场特征点。
本发明步骤1-2包括以下步骤:
步骤1-2-1,计算格点的涡度值:在全球风场格点数据中(全球数据规模为361×720,数据分辨率为0.5×0.5,其中赤道对应的行索引为180,本初子午线对应的列索引为0),格点涡度值计算公式如下所示:
Figure BDA0001645681850000043
其中,vorticityij为第i行第j列格点的涡度值,llon为当前纬度上相邻1°的经线之间距离,llat为相邻1°的纬线之间距离,vi,1为第i行第1列格点的经向风速,vi,719为第i行第719列格点的经向风速,vi,0为第i行第1列格点的经向风速,vi,718为第i行第718列格点的经向风速,vi,j+1为第i行第j+1列格点的经向风速,vi,j-1为第i行第j-1列格点的经向风速,ui-1,0为第i-1行第0列格点的纬向风速,ui+1,0为第i+1行第0列格点的纬向风速,ui-1,719为第i-1行第719列格点的纬向风速,ui+1,719为第i+1行第719列格点的纬向风速,ui-1,j为第i-1行第j列格点的纬向风速,ui+1,j为第i+1行第j列格点的纬向风速;
步骤1-2-2,设置涡度阈值,对风场特征点进行二次提取:涡度阈值的计算公式如下所示:
Figure BDA0001645681850000044
其中thresholdvorticity为涡度阈值,Φ为涡度阈值系数(根据多次提取结果,在最优结果中确定系数Φ的值,一般设置为3),sumvorticity为范围内格点涡度值之和,row为风场格点数据行数,rank为风场格点数据列数。
根据涡度阈值对初步提取的风场特征点进行二次提取,将涡度值大于涡度阈值的格点最终提取为风场特征点。
本发明步骤2包括以下步骤:
步骤2-1,槽线补充分析判别:在槽线处设定补充分析判别区域,根据该区域中风场特征点的数量,判断是否需对该条槽线进行补充分析;
步骤2-2,槽线订正判别:在槽线处设定订正区域,根据该区域内风场特征点在槽线周围的分布情况,判断是否需要订正该槽线;
本发明步骤2-1包括以下步骤:
步骤2-1-1,设定槽线补充分析判别区域:以槽线首尾端点之间连线为对角线,作槽线外接矩形A,将矩形A的对角线沿其两端点的方向各延长两倍,分别作矩形B和矩形C,令矩形A、B之间所夹区域为补充分析判别区域一Additional_Region_1,矩形A、C之间所夹区域为补充分析判别区域二Additional_Region_2;
步骤2-1-2,槽线补充分析判断:分别统计判别区域一和判别区域二内的风场特征点数量,若其中至少一个区域中风场特征点数量达到设定阈值thresholdadditional(一般设定为12),则判断需对该槽线进行补充分析,并将数量达到设定阈值的判别区域所对应的槽线端点,设置为槽线补充分析起点addStartpoint;
本发明步骤2-2包括以下步骤:
步骤2-2-1,设定槽线订正判别区域:遍历当前槽线上所有槽点,分别提取其中行索引最大点、行索引最小点以及列索引最大点、列索引最小的点(若存在多个槽点的索引均为最值,则任取其一;最值点可以重合,例如行索引最大点和列索引最大点可以为同一槽点);
将行索引最小点的行索引值减去2,得到订正判别区域的上边界点upPoint,将行索引最大点的行索引值加上2,得到下边界点downPoint,将列索引最小点的列索引值减去2,得到订正判别区域的左边界点leftPoint,将列索引最大点的列索引值加上2,得到右边界点rightPoint,根据上、下、左、右四个边界点构造该槽线的订正判别区域Correction_Region;
步骤2-2-2,槽线订正判断:建立特征点集F,存储订正判别区域Correction_Region中不在槽线上的风场特征点,并为此类特征点分别建立其属性表attributei(P1,l),其中i为此类风场特征点的索引,P1表示槽线上与该特征点相距最短的槽点的集合(若特征点存在多个距离最短槽点,则均需存储于集合P1中),l表示该最短距离;
将点集F中特征点按照最短距离l分类,在各类特征点中分别查找每两个特征点之间各自集合P1的交集,将该交集记为U,对集合U可能存在的情况进行分类讨论:
(1)若当前两个特征点的交集U为空集,则继续查找下两个特征点;
(2)若当前两个特征点的交集U中含有一个槽点,且两特征点关于该槽点对称,则在特征点集F中标记这两个特征点;
(3)若当前两个特征点的交集U中含有两个槽点,且该两特征点与此两槽点的连线互相垂直,则同样在特征点集F中标记这两个特征点(点集F中特征点可重复标记);
统计特征点集F中被标记的特征点数量,若大于F中特征点总数的一半,则无需对其进行订正,否则即判断该槽线需进行订正处理。
本发明步骤3包括以下步骤:
步骤3-1,补充分析点集提取:以补充分析判别区域对应的槽线端点作为补充分析起点,然后采用最小生成树聚类算法提取槽线补充分析所需的点集Add_Points;
步骤3-2,曲线拟合:根据槽线补充分析点集Add_Points内风场特征点的分布情况,采用四次多项式曲线拟合算法获得补充分析槽线;
本发明步骤3-1包括以下步骤:
步骤3-1-1:对于需进行补充分析的槽线,将槽线补充分析判别区域所对应的槽线端点作为槽线补充分析起点addStartpoint;
步骤3-1-2:根据风场特征点的坐标,计算两特征点之间的几何距离作为此两点之间的权值,以此构建加权连通图,用G=(Vs,Es)表示,其中Vs为风场特征点集,用于存储槽线的两个补充分析判别区域Additional_Region_1和Additional_Region_2重合部分之外的风场特征点,Es为特征点之间的加权边集,此外,设置空点集Vtree(初始为空集)和空边集Etree(初始为空集);
步骤3-1-3:将槽线补充分析起点addStartpoint加入点集Vtree中,在加权边集Es中遍历查找包含起点addStartpoint的权值最小边,并将最小边加入边集Etree中,同时将该边上不在Vtree中的特征点加入点集Vtree中;
步骤3-1-4:在加权边集Es中选取含有点集Vtree中特征点,且不包含于边集Etree的权值最小边emin(如存在多条符合条件的边,则任选其一),将该边加入边集Etree中,同时将边emin中不在点集Vtree中的特征点加入Vtree中;
步骤3-1-5:重复步骤3-1-4,直至点集Vtree=Vs
步骤3-1-6:设定聚类阈值cluster_thred=10,将边集Etree中权值大于cluster_thred的边删除,即将边集Etree分为多个子集,其中包含槽线补充分析起点addStartpoint的子集即为槽线补充分析点集Add_Points;
本发明步骤3-2包括以下步骤:
步骤3-2-1:设定拟合槽线多项式,如下所示:
y=a0+a1x+a2x2+a3x3+a4x4
其中x、y分别为拟合槽线上点的横坐标和纵坐标,a0、a1、a2、a3、a4分别为四次多项式的系数;
步骤3-2-2:计算在槽线补充分析点集Add_Points中,各风场特征点与拟合槽线的偏差平方和(设槽线补充分析点集中有n个风场特征点),计算公式如下所示:
Figure BDA0001645681850000071
其中,R2为偏差平方和,xi、yi分别表示索引为i的风场特征点的横坐标和纵坐标;
步骤3-2-3:将步骤3-2-2中公式的等式右边分别对系数a0、a1、a2、a3、a4求偏导,得到如下5个等式:
Figure BDA0001645681850000072
Figure BDA0001645681850000073
Figure BDA0001645681850000081
Figure BDA0001645681850000082
Figure BDA0001645681850000083
步骤3-2-4:对步骤3-2-3中5个等式的左边进行化简,并将其表示成矩阵形式,设槽线补充分析起点addStartpoint的坐标为(xas,yas),其中xas、yas分别为槽线补充分析起点的横坐标和纵坐标,在保证槽线补充分析起点位于拟合槽线的条件下,得到范德蒙德扩充矩阵方程,如下所示:
Figure BDA0001645681850000084
通过求解上述范德蒙德扩充矩阵方程矩阵,可得到拟合槽线上各槽点的横坐标、纵坐标,得到槽线补充分析结果。
本发明步骤4包括以下步骤:
步骤4-1:建立Laplacian变形框架:建立关于槽线的三角网格模型M=(V,E),其中V为槽线上槽点的集合,E为槽点之间边的集合,然后通过Laplacian矩阵L将槽线上各点几何坐标转换为Laplacian坐标;
步骤4-2:确定槽线订正锚点:建立槽点与风场特征点之间的对应关系,然后根据对应关系及权重具体确定槽线订正锚点的位置;
步骤4-3:求解优化方程:通过求解带位置约束的优化方程,获得订正之后各槽点的几何坐标,连接槽点并平滑完成槽线订正处理;
本发明步骤4-1包括以下步骤:
步骤4-1-1:对于单条槽线,建立关于槽线的三角网格模型M=(V,E),其中V为该槽线上槽点的集合,E为槽点之间边的集合;
步骤4-1-2:通过Laplacian转换矩阵L将槽线上各点的几何坐标转换为Laplacian坐标,计算公式如下所示:
LV=δ
其中,L为Laplacian转换矩阵,δ为经过转换后得到的槽点Laplacian坐标矩阵,而Laplacian转换矩阵L的形式如下:
Figure BDA0001645681850000091
其中,Lij为Laplacian转换矩阵L中第i行、第j列的元素,di为槽点集合V中第i个点的度,即该点邻接点的数量;
本发明步骤4-2包括以下步骤:
步骤4-2-1:提取特征点集F中未标记的风场特征点,并列举此类风场特征点的属性表attributei(P1,l)中点集P1所包含的全部槽点;
步骤4-2-2:为步骤4-2-1中列举的每个槽点分别建立点集P2,存储与该槽点相关的未标记特征点,反向映射槽点和风场特征点之间的关系;
步骤4-2-3:根据点集P2中包含特征点数量的情况,分情况建立槽点和风场特征点之间的对应关系:
若点集P2中仅含有一个风场特征点,则在该风场特征点与点集P2所对应的槽点之间建立对应关系;
若点集P2中含有两个以上风场特征点,则选取其中最短距离l最大的风场特征点(若存在多个风场特征点的l相同,则取其中涡度值vorticity最大的风场特征点),与点集P2所对应的槽点建立对应关系;
步骤4-2-4:在每组风场特征点和槽点的对应关系中,通过如下公式分别计算槽点和风场特征点在订正锚点确定过程中所占的比重:
Figure BDA0001645681850000092
其中,Weight(vorticity,αuv)为槽点和风场特征点在订正锚点确定过程中所占的比重,vorticity为涡度值,αuv分别为纬向风向偏转角度和经向风向偏转角度,vorticitysum、αu_sum、αv_sum分别为每组对应关系中,槽点和风场特征点的涡度值、纬向风向偏转角度之和、经向风向偏转角度之和;
步骤4-2-5:在每组对应关系中,根据槽点和风场特征点所占的比重Weight,确定槽线订正锚点的几何坐标(anchorPoint.X,anchorPoint.Y),如下所示:
Figure BDA0001645681850000101
Figure BDA0001645681850000102
其中,fp.X和fp.Y分别为风场特征点的横几何坐标、纵几何坐标,tp.X和tp.Y分别为槽点的横几何坐标、纵几何坐标,Weightfp、Weighttp分别为每组对应关系中风场特征点的比重和槽点的比重;
本发明步骤4-3包括以下步骤:
步骤4-3-1:加入槽线订正锚点后,将步骤4-1-2中的计算公式LV=δ变为如下所示:
Figure BDA0001645681850000103
其中,H为m×2阶矩阵,m为订正锚点个数,矩阵中每一行只有一个非零元素1,该元素表示订正锚点的权值,h为m×2阶矩阵,且hj=ωjUj,j=1,2,...,m,ωj即为权重值,hj即为矩阵h中的元素,Uj为槽线订正锚点的几何坐标,V'为n×2阶矩阵,其表示Laplacian变形后的槽点几何坐标,其第i行即表示Laplacian变形后的第i个槽点的几何坐标Vi'(xi,yi),i=1,2...,n;
步骤4-3-2:在步骤4-3-1中公式的等式两边同时左乘转置矩阵[L′]T,即可得到如下方程组:
Figure BDA0001645681850000111
步骤4-3-3:通过求逆得到步骤4-3-2中方程组的解,如下所示:
Figure BDA0001645681850000112
在求方程组解的过程中,首先将矩阵[L']TL'分解为上三角矩阵和下三角矩阵,然后结合位置约束条件对方程组进行多次迭代求解,最终获得较为精确的槽线订正处理后的槽点几何坐标,其中,位置约束条件如下所示:
Figure BDA0001645681850000113
步骤4-3-4:按照步骤4-3-1至步骤4-3-3,分别对槽点的横坐标、纵坐标进行求解,得到槽线订正处理后的槽点几何坐标,按序连接变形后的槽点并采用B样条曲线函数(参考文献:Wang W,Pottmann H,Liu Y.Fitting B-spline curves to point clouds bycurvature-based squared distance minimization[J].ACM Transactions on Graphics(ToG),2006,25(2):214-238)做平滑处理,得到槽线订正结果。
有益效果
本发明提出了一种基于风场数据的槽线修正方法,未来可根据该方法对已有的槽线自动分析方法进行改进和修正,建立更加准确的面向气象应用的槽线自动分析系统,减轻人工分析的压力,进一步提高天气图分析的实时性和准确性,完善天气图自动分析系统并推动其广泛使用。
附图说明
下面结合附图和具体实施方式对本发明做更进一步的具体说明,本发明的上述或其他方面的优点将会变得更加清楚。
图1为本发明流程图。
图2为角度坐标系示意图。
图3为槽线补充分析判别区域示意图。
图4为槽线订正判别区域示意图。
图5a为风场特征点相对槽线位置示意图(A、B两点)。
图5b为风场特征点相对槽线位置示意图(C、D两点)。
图6为补充分析点集提取示意图。
图7为Laplacian坐标示意图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
如图1所示,本发明提供了计算机气象软件中基于风场数据的槽线修正方法,具体包括如下步骤:
步骤1,提取风场特征点:获取风场格点数据,根据风场格点数据计算各格点经向和纬向上的风速,由风速计算结果得到各点风向角度,并映射到角度坐标系中,根据风向角度计算各点的纬向、经向风向偏转角度以及其涡度,然后分别设定经向、纬向风向偏转角度阈值和涡度阈值,在风场格点数据中提取风场特征点;
步骤2,判别槽线修正类型:在槽线周围设定槽线补充分析判别区域,根据该区域中风场特征点的数量,判断是否需进行槽线补充分析,如需要则执行步骤3;设定槽线订正判别区域,根据该区域中风场特征点与槽线的相对位置,判断是否需进行槽线订正处理,如需要则执行步骤4;
步骤3,槽线补充分析:根据步骤2所判别得到的槽线修正类型,对于需进行补充分析的槽线,首先以补充分析判别区域对应的槽线端点作为补充分析起点,然后在风场特征点中提取槽线补充分析所需的点集,最后由补充分析起点开始,在点集中通过曲线拟合完成槽线补充分析;
步骤4,槽线订正处理:根据步骤2所判别得到的槽线修正类型,对于需进行订正的槽线,首先存储网格信息并将槽线上槽点的几何坐标转换为Laplacian坐标,建立Laplacian坐标系,同时建立Laplacian坐标系,其中,Laplacian坐标的定义如下,
Figure BDA0001645681850000121
其中,δi即为顶点vi的Laplacian坐标,L()为Laplacian算子,vj为顶点vi的邻接点,N(i)为顶点vi的邻接点索引的集合,ωij为vj相对顶点vi的权值,目前常用的包括均匀权值、余切权值、正切权值等,本发明采用均匀权值,即ωij=1/di(di为vi点的度,即vi点邻接点的数量);
然后根据风场特征点和槽点确定槽线订正锚点,最后,通过求解带位置约束的优化方程,获得订正之后各槽点的几何坐标,连接槽点并平滑以完成槽线订正处理;
步骤5,根据步骤3、4的结果,输出槽线分析结果并显示在计算机屏幕上。
本发明步骤1包括以下步骤:
步骤1-1,初步提取风场特征点:于欧洲中期天气预报中心(ECMWF)官网下载获取风场格点数据,在风场格点数据中,根据其中各点的经向、纬向风速分量计算格点风向角度,并映射到角度坐标系中,然后对于全球范围内的离散风场格点数据,分别计算各格点的纬向、经向风向偏转角度,然后根据各格点的纬向、经向风向偏转角度,设置纬向、经向风向偏转角度阈值,对风场特征点进行初步提取;
步骤1-2,二次提取风场特征点:根据各格点的经向、纬向风速分量,计算格点的涡度值,然后根据格点涡度值设置涡度阈值,对风场特征点进行二次提取。
本发明步骤1-1包括以下步骤:
步骤1-1-1,通过如下公式计算格点风向角度:
Figure BDA0001645681850000131
其中,αij为第i行第j列的格点的风向角度,uij为该格点对应的纬向风速分量,vij为该格点对应的经向风速分量,由此分别计算各格点的风向角度;
步骤1-1-2,将格点风向角度映射至角度坐标系:将角度坐标系规定如下:定义正北为0°,正东为90°,正南为180°,正西为270°,如图2所示,根据经向、纬向的风速分量正负,将格点风向角度映射至角度坐标系中,公式如下所示:
Figure BDA0001645681850000132
映射完成之后,得到在角度坐标系下第i行第j列格点对应的风向角度值θij
步骤1-1-3,计算格点的纬向、经向风向偏转角度:在全球风场格点数据中(全球数据规模为361×720,数据分辨率为0.5×0.5,其中赤道对应的行索引为180,本初子午线对应的列索引为0),通过如下公式计算格点的纬向、经向风向偏转角度:
Figure BDA0001645681850000141
Figure BDA0001645681850000142
其中,αu_ij为第i行第j列格点的纬向风向偏转角度,αv_ij为第i行第j列格点的经向风向偏转角度,αi,1为第i行第1列格点的风向角度,αi,719为第i行第719列格点的风向角度,αi,j+1为第i行第j+1列格点的风向角度,αi+1,j为第i+1行第j列格点的风向角度;
步骤1-1-4,设置风向偏转角度阈值,对风场特征点进行初步提取:纬向、经向风向偏转角度阈值的计算公式如下所示:
Figure BDA0001645681850000143
Figure BDA0001645681850000144
其中,thresholdu、thresholdv分别为纬向、经向风向偏转角度阈值,
Figure BDA0001645681850000145
为风向偏转角度阈值系数,sumu和sumv分别为范围内格点纬向、经向风向偏转角度之和,row为风场格点数据行数,rank为风场格点数据列数;
规定若格点的纬向、经向风向偏转角度中,至少有一项大于其对应的风向偏转角度阈值,则将该格点初步提取为风场特征点。
本发明步骤1-2包括以下步骤:
步骤1-2-1,计算格点的涡度值:在全球风场格点数据中(全球数据规模为361×720,数据分辨率为0.5×0.5,其中赤道对应的行索引为180,本初子午线对应的列索引为0),格点涡度值计算公式如下所示:
Figure BDA0001645681850000151
其中,vorticityij为第i行第j列格点的涡度值,llon为当前纬度上相邻1°的经线之间距离,llat为相邻1°的纬线之间距离,vi,1为第i行第1列格点的经向风速,vi,719为第i行第719列格点的经向风速,vi,0为第i行第1列格点的经向风速,vi,718为第i行第718列格点的经向风速,vi,j+1为第i行第j+1列格点的经向风速,vi,j-1为第i行第j-1列格点的经向风速,ui-1,0为第i-1行第0列格点的纬向风速,ui+1,0为第i+1行第0列格点的纬向风速,ui-1,719为第i-1行第719列格点的纬向风速,ui+1,719为第i+1行第719列格点的纬向风速,ui-1,j为第i-1行第j列格点的纬向风速,ui+1,j为第i+1行第j列格点的纬向风速;
步骤1-2-2,设置涡度阈值,对风场特征点进行二次提取:涡度阈值的计算公式如下所示:
Figure BDA0001645681850000152
其中thresholdvorticity为涡度阈值,Φ为涡度阈值系数,sumvorticity为范围内格点涡度值之和,row为风场格点数据行数,rank为风场格点数据列数。
根据涡度阈值对初步提取的风场特征点进行二次提取,将涡度值大于涡度阈值的格点最终提取为风场特征点。
本发明步骤2包括以下步骤:
步骤2-1,槽线补充分析判别:在槽线处设定补充分析判别区域,根据该区域中风场特征点的数量,判断是否需对该条槽线进行补充分析;
步骤2-2,槽线订正判别:在槽线处设定订正区域,根据该区域内风场特征点在槽线周围的分布情况,判断是否需要订正该槽线;
本发明步骤2-1包括以下步骤:
步骤2-1-1,设定槽线补充分析判别区域:以槽线首尾端点之间连线为对角线,作槽线外接矩形A,将矩形A的对角线沿其两端点的方向各延长两倍,分别作矩形B和矩形C,令矩形A、B之间所夹区域为补充分析判别区域一Additional_Region_1,矩形A、C之间所夹区域为补充分析判别区域二Additional_Region_2,如图3所示;
步骤2-1-2,槽线补充分析判断:分别统计判别区域一和判别区域二内的风场特征点数量,若其中至少一个区域中风场特征点数量达到设定阈值thresholdadditional,则判断需对该槽线进行补充分析,并将数量达到设定阈值的判别区域所对应的槽线端点,设置为槽线补充分析起点addStartpoint;
本发明步骤2-2包括以下步骤:
步骤2-2-1,设定槽线订正判别区域:遍历当前槽线上所有槽点,分别提取其中行索引最大点、行索引最小点以及列索引最大点、列索引最小的点(若存在多个槽点的索引均为最值,则任取其一;最值点可以重合,例如行索引最大点和列索引最大点可以为同一槽点);
将行索引最小点的行索引值减去2,得到订正判别区域的上边界点upPoint,将行索引最大点的行索引值加上2,得到下边界点downPoint,将列索引最小点的列索引值减去2,得到订正判别区域的左边界点leftPoint,将列索引最大点的列索引值加上2,得到右边界点rightPoint,根据上、下、左、右四个边界点构造该槽线的订正判别区域Correction_Region,如图4所示;
步骤2-2-2,槽线订正判断:建立特征点集F,存储订正判别区域Correction_Region中不在槽线上的风场特征点,并为此类特征点分别建立其属性表attributei(P1,l),其中i为此类风场特征点的索引,P1表示槽线上与该特征点相距最短的槽点的集合(若特征点存在多个距离最短槽点,则均需存储于集合P1中),l表示该最短距离;
将点集F中特征点按照最短距离l分类,在各类特征点中分别查找每两个特征点之间各自集合P1的交集,将该交集记为U,对集合U可能存在的情况进行分类讨论:
(1)若当前两个特征点的交集U为空集,则继续查找下两个特征点;
(2)若当前两个特征点的交集U中含有一个槽点,且两特征点关于该槽点对称,如图5a中A、B两点所示,则在特征点集F中标记这两个特征点;
(3)若当前两个特征点的交集U中含有两个槽点,且该两特征点与此两槽点的连线互相垂直,如图5b中C、D两点所示,则同样在特征点集F中标记这两个特征点(点集F中特征点可重复标记);
统计特征点集F中被标记的特征点数量,若大于F中特征点总数的一半,则无需对其进行订正,否则即判断该槽线需进行订正处理。
本发明步骤3包括以下步骤:
步骤3-1,补充分析点集提取:以补充分析判别区域对应的槽线端点作为补充分析起点,然后采用最小生成树聚类算法提取槽线补充分析所需的点集Add_Points,如图6所示;
步骤3-2,曲线拟合:根据槽线补充分析点集Add_Points内风场特征点的分布情况,采用四次多项式曲线拟合算法获得补充分析槽线;
本发明步骤3-1包括以下步骤:
步骤3-1-1:对于需进行补充分析的槽线,将补充分析判别区域所对应的槽线端点作为补充分析起点addStartpoint;
步骤3-1-2:根据风场特征点的坐标,计算两特征点之间的几何距离作为此两点之间的权值,以此构建加权连通图,用G=(Vs,Es)表示,其中Vs为风场特征点集,用于存储槽线的两个补充分析判别区域Additional_Region_1和Additional_Region_2重合部分之外的风场特征点,Es为特征点之间的加权边集,此外,设置空点集Vtree(初始为空集)和空边集Etree(初始为空集);
步骤3-1-3:将槽线补充分析起点addStartpoint加入点集Vtree中,在加权边集Es中遍历查找包含起点addStartpoint的权值最小边,并将最小边加入边集Etree中,同时将该边上不在Vtree中的特征点加入点集Vtree中;
步骤3-1-4:在加权边集Es中选取含有点集Vtree中特征点,且不包含于边集Etree的权值最小边emin(如存在多条符合条件的边,则任选其一),将该边加入边集Etree中,同时将边emin中不在点集Vtree中的特征点加入Vtree中;
步骤3-1-5:重复步骤3-1-4,直至点集Vtree=Vs
步骤3-1-6:设定聚类阈值cluster_thred=10,将边集Etree中权值大于cluster_thred的边删除,即将边集Etree分为多个子集,其中包含补充分析起点addStartpoint的子集即为槽线补充分析点集Add_Points;
本发明步骤3-2包括以下步骤:
步骤3-2-1:设定拟合槽线多项式,如下所示:
y=a0+a1x+a2x2+a3x3+a4x4
其中x、y为拟合槽线上点的横、纵坐标,a0、a1、a2、a3、a4分别为四次多项式的系数;
步骤3-2-2:计算在槽线补充分析点集Add_Points中,各风场特征点与拟合槽线的偏差平方和(设槽线补充分析点集中有n个风场特征点),计算公式如下所示:
Figure BDA0001645681850000181
其中,R2为偏差平方和,xi、yi表示索引为i的风场特征点的横、纵坐标;
步骤3-2-3:将步骤3-2-2中的等式右边分别对系数a0、a1、a2、a3、a4求偏导,得到如下5个等式:
Figure BDA0001645681850000182
Figure BDA0001645681850000183
Figure BDA0001645681850000184
Figure BDA0001645681850000185
Figure BDA0001645681850000186
步骤3-2-4:对步骤3-2-3中5个等式的左边进行化简,并将其表示成矩阵形式,设槽线补充分析起点addStartpoint的坐标为(xas,yas),其中xas、yas分别为槽线补充分析起点的横、纵坐标,在保证槽线补充分析起点位于拟合槽线的条件下,得到范德蒙德扩充矩阵方程,如下所示:
Figure BDA0001645681850000191
通过求解上述范德蒙德扩充矩阵方程矩阵,可得到拟合槽线上各槽点的横、纵坐标,得到槽线补充分析结果。
本发明步骤4包括以下步骤:
步骤4-1:Laplacian变形框架建立:建立关于槽线的三角网格模型M=(V,E),其中V为槽线上槽点的集合,E为槽点之间边的集合,然后通过Laplacian矩阵L将槽线上各点几何坐标转换为Laplacian坐标;
步骤4-2:槽线订正锚点确定:建立槽点与风场特征点之间的对应关系,然后根据对应关系及权重具体确定槽线订正锚点的位置;
步骤4-3:优化方程求解:通过求解带位置约束的优化方程,获得订正之后各槽点的几何坐标,连接槽点并平滑完成槽线订正处理;
本发明步骤4-1包括以下步骤:
步骤4-1-1:对于单条槽线,建立关于槽线的三角网格模型M=(V,E),其中V为该槽线上槽点的集合,E为槽点之间边的集合;
步骤4-1-2:通过Laplacian转换矩阵L将槽线上各点的几何坐标转换为Laplacian坐标,Laplacian坐标示意图如图7所示,计算公式如下所示:
LV=δ
其中,L为Laplacian转换矩阵,δ为经过转换后得到的槽点Laplacian坐标矩阵,而Laplacian转换矩阵L的形式如下:
Figure BDA0001645681850000192
其中,Lij为Laplacian转换矩阵L中第i行、第j列的元素,di为槽点集合V中第i个点的度,即该点邻接点的数量;
本发明步骤4-2包括以下步骤:
步骤4-2-1:提取特征点集F中未标记的风场特征点,并列举此类风场特征点的属性表attributei(P1,l)中点集P1所包含的全部槽点;
步骤4-2-2:为步骤4-2-1中列举的每个槽点分别建立点集P2,存储与该槽点相关的未标记特征点,反向映射槽点和风场特征点之间的关系;
步骤4-2-3:根据点集P2中包含特征点数量的情况,分情况建立槽点和风场特征点之间的对应关系:
若点集P2中仅含有一个风场特征点,则在该风场特征点与点集P2所对应的槽点之间建立对应关系;
若点集P2中含有多个风场特征点,则选取其中最短距离l最大的风场特征点(若存在多个风场特征点的l相同,则取其中涡度值vorticity最大的风场特征点),与点集P2所对应的槽点建立对应关系;
步骤4-2-4:在每组风场特征点和槽点的对应关系中,通过如下公式分别计算槽点和风场特征点在订正锚点确定过程中所占的比重:
Figure BDA0001645681850000201
其中,Weight(vorticity,αuv)为槽点和风场特征点在订正锚点确定过程中所占的比重,vorticity为涡度值,αuv为纬向,经向风向偏转角度,vorticitysum、αu_sum、αv_sum分别为每组对应关系中,槽点和风场特征点的涡度值、纬向风向偏转角度、经向风向偏转角度之和;
步骤4-2-5:在每组对应关系中,根据槽点和风场特征点所占的比重Weight,确定槽线订正锚点的几何坐标(anchorPoint.X,anchorPoint.Y),如下所示:
Figure BDA0001645681850000202
Figure BDA0001645681850000211
其中,fp.X和fp.Y分别为风场特征点的横、纵几何坐标,tp.X和tp.Y分别为槽点的横、纵几何坐标,Weightfp、Weighttp分别为每组对应关系中风场特征点和槽点的比重;
本发明步骤4-3包括以下步骤:
步骤4-3-1:加入槽线订正锚点后,将步骤4-1-2中的计算公式LV=δ变为如下所示:
Figure BDA0001645681850000212
其中,H为m×2阶矩阵,m为订正锚点个数,矩阵中每一行只有一个非零元素1,该元素表示订正锚点的权值,h为m×2阶矩阵,且hj=ωjUj,j=1,2,...,m,ωj即为权重值,hj即为矩阵h中的元素,Uj为槽线订正锚点的几何坐标,V'为n×2阶矩阵,其表示Laplacian变形后的槽点几何坐标,其第i行即表示Laplacian变形后的第i个槽点的几何坐标Vi'(xi,yi),i=1,2...,n;
步骤4-3-2:在步骤4-3-1的等式两边同时左乘转置矩阵[L′]T,即可得到如下方程组:
Figure BDA0001645681850000213
步骤4-3-3:通过求逆得到步骤4-3-2中方程组的解,如下所示:
Figure BDA0001645681850000214
在求方程组解的过程中,首先将矩阵[L']TL'分解为上三角矩阵和下三角矩阵,然后结合位置约束条件对方程组进行多次迭代求解,最终获得较为精确的槽线订正处理后的槽点几何坐标,其中,位置约束条件如下所示:
Figure BDA0001645681850000221
步骤4-3-4:按照步骤4-3-1至步骤4-3-3,分别对槽点的横、纵坐标进行求解,得到槽线订正处理后的槽点几何坐标,按序连接变形后的槽点并采用B样条曲线函数(参考文献:Wang W,Pottmann H,Liu Y.Fitting B-spline curves to point clouds bycurvature-based squared distance minimization[J].ACM Transactions on Graphics(ToG),2006,25(2):214-238)做平滑处理,得到槽线订正结果。
本发明提供了一种计算机气象软件中基于风场数据的槽线修正方法,具体实现该技术方案的方法和途径很多,以上所述仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (7)

1.一种计算机气象软件中基于风场数据的槽线修正方法,其特征在于,包括以下步骤:
步骤1,提取风场特征点:获取风场格点数据,分别设定经向风向偏转角度阈值、纬向风向偏转角度阈值和涡度阈值,在风场格点数据中提取风场特征点;
步骤2,判别槽线修正类型:判断是否需进行槽线补充分析,如需要则执行步骤3;然后判断是否需进行槽线订正处理,如需要则执行步骤4;
步骤3,槽线补充分析:在风场特征点中提取槽线补充分析所需的点集,在点集中通过曲线拟合完成槽线补充分析;
步骤4,槽线订正处理:存储风场格点数据中所包含的网格信息,网格信息即各格点几何坐标及其之间的距离关系,并将槽线上槽点的几何坐标转换为Laplacian坐标,同时建立Laplacian坐标系,其中,Laplacian坐标的定义如下,
Figure FDA0003496351330000011
其中,δi即为顶点vi的Laplacian坐标,L()为Laplacian算子,vj为顶点vi的邻接点,N(i)为顶点vi的邻接点索引的集合,ωij为vj相对顶点vi的权值;
根据风场特征点和槽点确定槽线订正锚点,求解订正之后各槽点的几何坐标,连接槽点并平滑以完成槽线订正处理;
步骤5,根据步骤3和步骤4的结果,输出槽线分析结果并显示在计算机屏幕上;
步骤1包括以下步骤:
步骤1-1,初步提取风场特征点:获取风场格点数据,根据其中各点的经向、纬向风速分量计算格点风向角度,并映射到角度坐标系中,对于全球范围内的离散风场格点数据,分别计算各格点的纬向、经向风向偏转角度,然后根据各格点的纬向、经向风向偏转角度,设置纬向、经向风向偏转角度阈值,对风场特征点进行初步提取;
步骤1-2,二次提取风场特征点:根据各格点的经向、纬向风速分量,计算格点的涡度值,根据格点涡度值设置涡度阈值,对风场特征点进行二次提取;
步骤1-1包括以下步骤:
步骤1-1-1,通过如下公式计算格点风向角度:
Figure FDA0003496351330000012
其中,αij为第i行第j列的格点的风向角度,uij为该格点对应的纬向风速分量,vij为该格点对应的经向风速分量;
步骤1-1-2,将格点风向角度映射至角度坐标系:将角度坐标系规定如下:定义正北为0°,正东为90°,正南为180°,正西为270°,根据经向、纬向的风速分量正负,将格点风向角度映射至角度坐标系中,公式如下所示:
Figure FDA0003496351330000021
映射完成之后,得到在角度坐标系下第i行第j列格点对应的风向角度值θij
步骤1-1-3,计算格点的纬向、经向风向偏转角度:在全球风场格点数据中,通过如下公式计算格点的纬向风向偏转角度和经向风向偏转角度:
Figure FDA0003496351330000022
Figure FDA0003496351330000023
其中,全球风场格点数据规模为361×720,数据分辨率为0.5×0.5,其中赤道对应的行索引为180,本初子午线对应的列索引为0,αu_ij为第i行第j列格点的纬向风向偏转角度,αv_ij为第i行第j列格点的经向风向偏转角度,αi,1为第i行第1列格点的风向角度,αi,719为第i行第719列格点的风向角度,αi,j+1为第i行第j+1列格点的风向角度,αi+1,j为第i+1行第j列格点的风向角度;
步骤1-1-4,设置风向偏转角度阈值,对风场特征点进行初步提取:纬向风向偏转角度阈值、经向风向偏转角度阈值的计算公式如下所示:
Figure FDA0003496351330000024
Figure FDA0003496351330000031
其中,thresholdu、thresholdv分别为纬向风向偏转角度阈值和经向风向偏转角度阈值,
Figure FDA0003496351330000032
为风向偏转角度阈值系数,sumu和sumv分别为范围内格点纬向风向偏转角度之和、经向风向偏转角度之和,row为风场格点数据行数,rank为风场格点数据列数;
若格点的纬向、经向风向偏转角度中,至少有一项大于其对应的风向偏转角度阈值,则将该格点初步提取为风场特征点;
步骤1-2包括以下步骤:
步骤1-2-1,计算格点的涡度值:在全球风场格点数据中,格点涡度值计算公式如下所示:
Figure FDA0003496351330000033
其中,vorticityij为第i行第j列格点的涡度值,llon为当前纬度上相邻1°的经线之间距离,llat为相邻1°的纬线之间距离,vi,1为第i行第1列格点的经向风速,vi,719为第i行第719列格点的经向风速,vi,0为第i行第1列格点的经向风速,vi,718为第i行第718列格点的经向风速,vi,j+1为第i行第j+1列格点的经向风速,vi,j-1为第i行第j-1列格点的经向风速,ui-1,0为第i-1行第0列格点的纬向风速,ui+1,0为第i+1行第0列格点的纬向风速,ui-1,719为第i-1行第719列格点的纬向风速,ui+1,719为第i+1行第719列格点的纬向风速,ui-1,j为第i-1行第j列格点的纬向风速,ui+1,j为第i+1行第j列格点的纬向风速;
步骤1-2-2,设置涡度阈值,对风场特征点进行二次提取:涡度阈值的计算公式如下所示:
Figure FDA0003496351330000034
其中thresholdvorticity为涡度阈值,Φ为涡度阈值系数,sumvorticity为范围内格点涡度值之和;
根据涡度阈值对初步提取的风场特征点进行二次提取,将涡度值大于涡度阈值的格点最终提取为风场特征点。
2.根据权利要求1所述的方法,其特征在于,步骤2包括以下步骤:
步骤2-1,槽线补充分析判别:在槽线处设定槽线补充分析判别区域,根据该区域中风场特征点的数量,判断是否需对该条槽线进行补充分析;
步骤2-2,槽线订正判别:在槽线处设定槽线订正判别区域,根据该区域内风场特征点在槽线周围的分布情况,判断是否需要订正该槽线。
3.根据权利要求2所述的方法,其特征在于,步骤2-1包括以下步骤:
步骤2-1-1,设定槽线补充分析判别区域:以槽线首尾端点之间连线为对角线,作槽线外接矩形A,将矩形A的对角线沿其两端点的方向各延长两倍,分别作矩形B和矩形C,令矩形A、B之间所夹区域为补充分析判别区域一Additional_Region_1,矩形A、C之间所夹区域为补充分析判别区域二Additional_Region_2;
步骤2-1-2,槽线补充分析判断:分别统计判别区域一和判别区域二内的风场特征点数量,若其中至少一个区域中风场特征点数量达到设定阈值thresholdadditional,则判断需对该槽线进行补充分析,并将数量达到设定阈值的判别区域所对应的槽线端点,设置为槽线补充分析起点addStartpoint。
4.根据权利要求3所述的方法,其特征在于,步骤2-2包括以下步骤:
步骤2-2-1,设定槽线订正判别区域:遍历当前槽线上所有槽点,分别提取其中行索引最大点、行索引最小点以及列索引最大点、列索引最小的点;
将行索引最小点的行索引值减去2,得到订正判别区域的上边界点upPoint,将行索引最大点的行索引值加上2,得到下边界点downPoint,将列索引最小点的列索引值减去2,得到订正判别区域的左边界点leftPoint,将列索引最大点的列索引值加上2,得到右边界点rightPoint,根据上、下、左、右四个边界点构造该槽线的订正判别区域Correction_Region;
步骤2-2-2,槽线订正判断:建立特征点集F,存储订正判别区域Correction_Region中不在槽线上的风场特征点,并为此类特征点分别建立其属性表attributei(P1,l),其中i为此类风场特征点的索引,P1表示槽线上与该特征点相距最短的槽点集合,l表示该最短距离;
将点集F中特征点按照最短距离l分类,在各类特征点中分别查找每两个特征点之间各自集合P1的交集,将该交集记为U,对集合U可能存在的情况进行如下分类讨论:
若当前两个特征点的交集U为空集,则继续查找下两个特征点;
若当前两个特征点的交集U中含有一个槽点,且两特征点关于该槽点对称,则在特征点集F中标记这两个特征点;
若当前两个特征点的交集U中含有两个槽点,且该两特征点与此两槽点的连线互相垂直,则同样在特征点集F中标记这两个特征点;
统计特征点集F中被标记的特征点数量,若大于F中特征点总数的一半,则不进行订正,否则判断该槽线需进行订正处理。
5.根据权利要求4所述的方法,其特征在于,步骤3包括以下步骤:
步骤3-1,补充分析点集提取:以槽线补充分析判别区域对应的槽线端点作为槽线补充分析起点,采用最小生成树聚类算法提取槽线补充分析所需的点集Add_Points;
步骤3-2,曲线拟合:根据槽线补充分析点集Add_Points内风场特征点的分布情况,采用四次多项式曲线拟合算法获得补充分析槽线;
其中,步骤3-1包括以下步骤:
步骤3-1-1:对于需进行补充分析的槽线,将槽线补充分析判别区域所对应的槽线端点作为槽线补充分析起点addStartpoint;
步骤3-1-2:根据风场特征点的坐标,计算两特征点之间的几何距离作为此两点之间的权值,以此构建加权连通图,用G=(Vs,Es)表示,其中Vs为风场特征点集,用于存储槽线的两个补充分析判别区域Additional_Region_1和Additional_Region_2重合部分之外的风场特征点,Es为特征点之间的加权边集,设置初始为空的空点集Vtree和空边集Etree
步骤3-1-3:将槽线补充分析起点addStartpoint加入点集Vtree中,在加权边集Es中遍历查找包含起点addStartpoint的权值最小边,并将最小边加入边集Etree中,同时将该边上不在Vtree中的特征点加入点集Vtree中;
步骤3-1-4:在加权边集Es中选取含有点集Vtree中特征点,且不包含于边集Etree的权值最小边emin,将该边加入边集Etree中,同时将边emin中不在点集Vtree中的特征点加入Vtree中;
步骤3-1-5:重复步骤3-1-4,直至点集Vtree=Vs
步骤3-1-6:设定聚类阈值cluster_thred=10,将边集Etree中权值大于cluster_thred的边删除,即将边集Etree分为多个子集,其中包含槽线补充分析起点addStartpoint的子集即为槽线补充分析点集Add_Points;
其中,步骤3-2包括以下步骤:
步骤3-2-1:设定拟合槽线多项式,如下所示:
y=a0+a1x+a2x2+a3x3+a4x4
其中x、y分别为拟合槽线上点的横坐标和纵坐标,a0、a1、a2、a3、a4分别为四次多项式的系数;
步骤3-2-2:设槽线补充分析点集中有n个风场特征点,计算在槽线补充分析点集Add_Points中,各风场特征点与拟合槽线的偏差平方和,计算公式如下所示:
Figure FDA0003496351330000061
其中,R2为偏差平方和,xi、yi分别表示索引为i的风场特征点的横坐标和纵坐标;
步骤3-2-3:将步骤3-2-2中公式的等式右边分别对系数a0、a1、a2、a3、a4求偏导,得到如下5个等式:
Figure FDA0003496351330000062
Figure FDA0003496351330000063
Figure FDA0003496351330000064
Figure FDA0003496351330000065
Figure FDA0003496351330000071
步骤3-2-4:对步骤3-2-3中5个等式的左边进行化简,并将其表示成矩阵形式,设槽线补充分析起点addStartpoint的坐标为(xas,yas),其中xas、yas分别为槽线补充分析起点的横坐标和纵坐标,在保证槽线补充分析起点位于拟合槽线的条件下,得到范德蒙德扩充矩阵方程,如下所示:
Figure FDA0003496351330000072
通过求解上述范德蒙德扩充矩阵方程矩阵,得到拟合槽线上各槽点的横坐标、纵坐标,得到槽线补充分析结果。
6.根据权利要求5所述的方法,其特征在于,步骤4包括以下步骤:
步骤4-1:建立Laplacian变形框架:建立关于槽线的三角网格模型M=(V,E),其中V为槽线上槽点的集合,E为槽点之间边的集合,通过Laplacian矩阵L将槽线上各点几何坐标转换为Laplacian坐标;
步骤4-2:确定槽线订正锚点:建立槽点与风场特征点之间的对应关系,根据对应关系及权重具体确定槽线订正锚点的位置;
步骤4-3:求解优化方程:通过求解带位置约束的优化方程,获得订正之后各槽点的几何坐标,连接槽点并平滑完成槽线订正处理。
7.根据权利要求6所述的方法,其特征在于,步骤4-1包括以下步骤:
步骤4-1-1:对于单条槽线,建立关于槽线的三角网格模型M=(V,E),其中V为该槽线上槽点的集合,E为槽点之间边的集合;
步骤4-1-2:通过Laplacian转换矩阵L将槽线上各点的几何坐标转换为Laplacian坐标,计算公式如下所示:
LV=δ,
其中,L为Laplacian转换矩阵,δ为经过转换后得到的槽点Laplacian坐标矩阵,而Laplacian转换矩阵L的形式如下:
Figure FDA0003496351330000081
其中,Lij为Laplacian转换矩阵L中第i行、第j列的元素,di为槽点集合V中第i个点的度,即该点邻接点的数量;
其中,步骤4-2包括以下步骤:
步骤4-2-1:提取特征点集F中未标记的风场特征点,并列举此类风场特征点的属性表attributei(P1,l)中点集P1所包含的全部槽点;
步骤4-2-2:为步骤4-2-1中列举的每个槽点分别建立点集P2,存储与该槽点相关的未标记特征点,反向映射槽点和风场特征点之间的关系;
步骤4-2-3:根据点集P2中包含特征点数量的情况,分情况建立槽点和风场特征点之间的对应关系:
若点集P2中仅含有一个风场特征点,则在该风场特征点与点集P2所对应的槽点之间建立对应关系;
若点集P2中含有两个以上风场特征点,则选取其中最短距离l最大的风场特征点,与点集P2所对应的槽点建立对应关系;
步骤4-2-4:在每组风场特征点和槽点的对应关系中,通过如下公式分别计算槽点和风场特征点在订正锚点确定过程中所占的比重:
Figure FDA0003496351330000082
其中,Weight(vorticity,αuv)为槽点和风场特征点在订正锚点确定过程中所占的比重,vorticity为涡度值,αuv分别为纬向风向偏转角度和经向风向偏转角度,vorticitysum、αu_sum、αv_sum分别为每组对应关系中,槽点和风场特征点的涡度值、纬向风向偏转角度之和、经向风向偏转角度之和;
步骤4-2-5:在每组对应关系中,根据槽点和风场特征点所占的比重Weight,确定槽线订正锚点的几何坐标(anchorPoint.X,anchorPoint.Y),如下所示:
Figure FDA0003496351330000091
Figure FDA0003496351330000092
其中,fp.X和fp.Y分别为风场特征点的横几何坐标、纵几何坐标,tp.X和tp.Y分别为槽点的横几何坐标、纵几何坐标,Weightfp、Weighttp分别为每组对应关系中风场特征点的比重和槽点的比重;
其中,步骤4-3包括以下步骤:
步骤4-3-1:加入槽线订正锚点后,将步骤4-1-2中的计算公式LV=δ变为如下所示:
Figure FDA0003496351330000093
其中,H为m×2阶矩阵,m为订正锚点个数,矩阵中每一行只有一个非零元素1,该元素表示订正锚点的权值,h为m×2阶矩阵,且hj=ωjUj,j=1,2,...,m,ωj即为权重值,hj即为矩阵h中的元素,Uj为槽线订正锚点的几何坐标,V'为n×2阶矩阵,其表示Laplacian变形后的槽点几何坐标,其第i行即表示Laplacian变形后的第i个槽点的几何坐标Vi'(xi,yi),i=1,2...,n;
步骤4-3-2:在步骤4-3-1中公式的等式两边同时左乘转置矩阵[′]T,得到如下方程组:
Figure FDA0003496351330000101
步骤4-3-3:通过求逆得到步骤4-3-2中方程组的解,如下所示:
Figure FDA0003496351330000102
在求方程组解的过程中,将矩阵[L']TL'分解为上三角矩阵和下三角矩阵,然后结合位置约束条件对方程组进行多次迭代求解,最终获得槽线订正处理后的槽点几何坐标,其中,位置约束条件如下所示:
Figure FDA0003496351330000103
步骤4-3-4:按照步骤4-3-1至步骤4-3-3,分别对槽点的横坐标、纵坐标进行求解,得到槽线订正处理后的槽点几何坐标,按序连接变形后的槽点并采用B样条曲线函数做平滑处理,得到槽线订正结果。
CN201810400897.5A 2018-04-28 2018-04-28 一种计算机气象软件中基于风场数据的槽线修正方法 Active CN108875127B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810400897.5A CN108875127B (zh) 2018-04-28 2018-04-28 一种计算机气象软件中基于风场数据的槽线修正方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810400897.5A CN108875127B (zh) 2018-04-28 2018-04-28 一种计算机气象软件中基于风场数据的槽线修正方法

Publications (2)

Publication Number Publication Date
CN108875127A CN108875127A (zh) 2018-11-23
CN108875127B true CN108875127B (zh) 2022-04-19

Family

ID=64326891

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810400897.5A Active CN108875127B (zh) 2018-04-28 2018-04-28 一种计算机气象软件中基于风场数据的槽线修正方法

Country Status (1)

Country Link
CN (1) CN108875127B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110346518B (zh) * 2019-07-25 2021-06-15 中南大学 一种交通排放污染可视化预警方法及其系统
CN110346517B (zh) * 2019-07-25 2021-06-08 中南大学 一种智慧城市工业大气污染可视化预警方法及其系统
CN112765832B (zh) * 2021-02-02 2022-05-06 南京信息工程大学 一种欧亚大陆冷锋自动识别订正方法
CN113486844B (zh) * 2021-07-26 2022-04-15 中科三清科技有限公司 水平风切变位置判断方法、装置、电子设备及存储介质

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7075493B2 (en) * 2002-05-01 2006-07-11 The Regents Of The University Of Michigan Slot antenna
KR101538933B1 (ko) * 2013-11-07 2015-07-23 재단법인 한국형수치예보모델개발사업단 수치일기예보모델의 데이터를 육면 격자 프레임에 가시화하는 방법 및 이를 수행하는 하드웨어 장치
CN104898186B (zh) * 2015-06-09 2017-05-03 天津大学 一种用于槽脊线的特征点提取及自动绘制的方法
CN104951624B (zh) * 2015-07-13 2017-10-17 中国人民解放军理工大学 一种计算机气象软件中基于风场数据的槽线自动绘制方法
US9842421B2 (en) * 2015-07-23 2017-12-12 Pixar Method and system for vorticle fluid simulation
CN106909788B (zh) * 2017-02-28 2019-02-22 中国人民解放军理工大学 计算机气象软件中基于位势高度数据的槽线自动绘制方法
CN106950612B (zh) * 2017-03-14 2019-07-23 天津大学 一种用于气象学中自动识别并绘制冷锋的方法

Also Published As

Publication number Publication date
CN108875127A (zh) 2018-11-23

Similar Documents

Publication Publication Date Title
CN108875127B (zh) 一种计算机气象软件中基于风场数据的槽线修正方法
CN111640125B (zh) 基于Mask R-CNN的航拍图建筑物检测和分割方法及装置
CN103727930B (zh) 一种基于边缘匹配的激光测距仪与相机相对位姿标定方法
JP7048225B2 (ja) 建物領域抽出用の学習済みモデル
CN107153822A (zh) 一种基于深度学习的半自动图像精标注方法
CN102592268A (zh) 一种分割前景图像的方法
CN102938066A (zh) 一种基于多元数据重建建筑物外轮廓多边形的方法
CN103077555B (zh) 一种三维模型构成的自动标注方法
CN102324102A (zh) 一种图像场景空洞区域结构和纹理信息自动填补方法
CN105574527A (zh) 一种基于局部特征学习的快速物体检测方法
CN104134220A (zh) 一种像方空间一致性的低空遥感影像高精度匹配方法
CN113689445B (zh) 结合语义分割与边缘检测的高分辨率遥感建筑物提取方法
CN104951624A (zh) 一种计算机气象软件中基于风场数据的槽线自动绘制方法
CN106485676A (zh) 一种基于稀疏编码的LiDAR点云数据修复方法
CN101930483A (zh) 应用参数化设计模型化简数字地图居民地多边形的方法
CN101833790B (zh) 一种基于波动方程的各向异性的四边形网格生成方法
CN106340010A (zh) 一种基于二阶轮廓差分的角点检测方法
CN106780450A (zh) 一种基于低秩多尺度融合的图像显著性检测方法
CN106599810A (zh) 一种基于栈式自编码的头部姿态估计方法
CN104036550A (zh) 基于形状语义的建筑立面激光雷达点云解译与重建的方法
CN103839274B (zh) 一种基于几何比例关系的扩展目标跟踪方法
CN105107901B (zh) 船舶双向曲率板自动成型加工路径确定及成型效果检测方法
CN104457691B (zh) 一种建筑物主体高程信息获取方法
CN102982567A (zh) 一种基于统计分析的变形体碰撞检测剔除方法
CN102663815B (zh) 一种基于水平集的lod2建筑物模型构建方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant