[发明专利]计算机气象软件中基于位势高度数据的槽线自动绘制方法

摘要

本发明公开了计算机气象软件中基于位势高度数据的槽线自动绘制方法，可以实现天气图中气象槽线的自动提取绘制，包括以下步骤：曲率值计算：计算各网格点经向和纬向上的地转风速，并由此计算得到网格点处地转风的风向角度，在风速方向上对地转风风向角度求导得到曲率值。槽点追踪：提取局部曲率极小值作为槽线追踪起点，由槽线追踪起点开始，根据各网格点的曲率值较小邻接点的数量等约束条件对槽点进行追踪和连接，得到原始槽线。原始槽线处理：对提取出的原始槽线进行噪声剔除，平滑处理，生成满足气象业务需求的槽线。

主权项

1.计算机气象软件中基于位势高度数据的槽线自动绘制方法，其特征在于，包括以下步骤：步骤1，计算曲率值：根据位势高度网格数据计算各网格点纬向和经向上的地转风风速，由风速计算结果得到地转风风向角度，并映射到角度坐标系中，在网格点风速方向上对地转风风向角度求导得到网格点的曲率值，该曲率值视为过该网格点的位势高度等值线在该点处的曲率值，最后计算各网格点的曲率值较小邻接点的数量；步骤2，槽点追踪：提取局部曲率极小值所在的网格点作为槽线追踪起点，从追踪起点开始，根据各网格点的曲率值较小邻接点的数量和相应约束条件，先沿曲率值较大的方向进行槽点的追踪和连接，当该方向追踪结束后，由槽线追踪起点开始对另一个曲率值较大的方向进行槽点的追踪和连接，得到原始槽线；步骤3，处理原始槽线：对步骤2追踪得到的原始槽线，进行噪声剔除和平滑处理，生成满足气象业务需求的槽线；步骤4，根据步骤3的结果，自动提取绘制位势高度场的槽线并显示在计算机屏幕上；步骤1包括以下步骤：步骤1‑1，计算纬向和经向上的地转风风速：在位势高度场网格数据中，在经向即垂直方向上对位势高度值求偏导，计算得到纬向上的地转风风速；在纬向即水平方向上对位势高度值求偏导，计算得到经向上的地转风风速；步骤1‑2，计算地转风风向角度：遍历位势高度场网格数据，根据纬向和经向上的地转风风速计算得到各网格点的地转风风向角度，映射到角度坐标系中；步骤1‑3，计算网格点的曲率值：利用双线性插值法得到网格点地转风风向上L/2处一网格点的风向，其中L为网格点距，在网格点的地转风方向上对地转风风向角度求偏导得到该网格点的曲率值；步骤1‑4，计算得到各网格点的曲率值较小邻接点的数量：遍历位势高度场网格数据中的网格点，比较该点与其周围八个邻接点的曲率值，得到每个网格点的曲率值较小邻接点数量；步骤1‑1包括以下步骤：步骤1‑1‑1，通过如下公式计算科里奥利参数：其中，f为该网格点对应的科里奥利参数，ω为地球自传角速度，为地理纬度，根据各网格点的地理纬度不同，分别计算出其所对应的科里奥利参数；步骤1‑1‑2，计算纬向的地转风风速：根据各网格点的科里奥利参数f，在经向即垂直方向上对位势高度值求偏导，乘上科里奥利参数f的‑1次方，对计算结果取相反数，得到纬向上的地转风风速，其计算公式如下所示：其中，u_g为纬向地转风风速，φ为位势高度值，y表示垂直方向，表示位势高度在垂直方向上的偏导数值；对于位势高度场网格数据中非首行和非末行的数据网格点，即行标大于1且小于数据行数Row的数据点，采用中央差分法计算位势高度在垂直方向上的偏导数值计算公式如下所示：其中，φi，j为第i行第j列的数据网格点的位势高度值，φi‑1，j为第i‑1行第j列的数据网格点的位势高度值，φi+1，j为第i+1行第j列的数据网格点的位势高度值，y表示垂直方向；对于位势高度场网格数据中的首行和末行数据，分别采用后项差分法和前项差分法计算位势高度在垂直方向上的偏导数值计算公式如下所示：其中，φRow，j为第Row行第j列的数据网格点的位势高度值，φRow‑1，j为第Row‑1行第j列的数据网格点的位势高度值；在作差分计算时，用行标较小网格点的位势高度值减去行标较大网格点的位势高度值；将各网格点位势高度在垂直方向上的偏导数值与科里奥利参数的‑1次方相乘，并将计算结果取相反数，得到纬向地转风风速；步骤1‑1‑3，计算经向的地转风风速：根据各网格点的科里奥利参数f，在纬向即水平方向上对位势高度值求偏导，乘上科里奥利参数f的‑1次方，得到经向上的地转风风速，其计算公式如下所示：其中，v_g为经向地转风风速，f为科里奥利参数，x表示水平方向，表示位势高度在水平方向上的偏导数值；对于位势高度场网格数据中非首列和非末列的数据网格点，即列标大于1且小于数据列数Rank的数据点，采用中央差分法计算位势高度在水平方向上的偏导数值计算公式如下所示：其中，φi，j+1为第i行第j+1列的数据网格点的位势高度值，φi，j‑1为第i行第j‑1列的数据网格点的位势高度值；对于网格数据中的首列和末列数据，分别采用后项差分法和前项差分法计算位势高度在水平方向上的偏导数值计算公式如下所示：φi，Rank为第i行第Rank列的数据网格点的位势高度值；在作差分计算时，用列标较大网格点的位势高度值减去列标较小网格点的位势高度值；将各网格点位势高度在水平方向上的偏导数值与科里奥利参数的‑1次方相乘，得到经向地转风风速；步骤1‑2包括以下步骤：步骤1‑2‑1，根据经向及纬向的地转风风速计算地转风风向角度：计算公式如下所示：其中，θi,j表示i行j列上的地转风风向角度，ugi,j为i行j列上的纬向地转风风速，vgi，j为i行j列上的经向地转风风速；步骤1‑2‑2，将地转风风向角度映射至角度坐标系：将角度坐标系规定如下：定义正北为0°，正东为90°，正南为180°，正西为270°，根据经向、纬向的地转风风速的正负，将地转风风向角度映射至角度坐标系中，公式如下所示：映射完成之后，得到在统一角度坐标系下全部数据网格点的地转风风向角度值θi,j；步骤1‑3包括以下步骤：步骤1‑3‑1，计算非边界网格点的曲率值：利用双线性插值法得到网格点地转风风向上一定距离处两点处的地转风风向角度，在网格点地转风风向上进行中央差分运算，得到非边界网格点处的曲率值；步骤1‑3‑2，计算边界网格点的曲率值：利用双线性插值法得到网格点地转风风向上一定距离处一点处的地转风风向角度，在网格点地转风风向上进行前项或后项差分运算，得到边界网格点处的曲率值；步骤1‑3‑1包括以下步骤：步骤1‑3‑1‑1，插值计算邻近两点的地转风风向角度：利用双线性插值法进行计算，网格点距为L，网格点距即相邻数据网格点之间的距离，其中，A，D，H，F，I，J，K为网格点，A为待求曲率值的网格点，D、H、F、K、I、J分别为A点的右侧网格点、右上方网格点、上方网格点、左侧网格点、左下方网格点和下方网格点，E、G、B分别为AD上的点、FH上的点和EG上的点，AB和AC的长度均为α为A点的地转风风向角度，E点的纬向地转风风速和经向地转风速由A、D两点的数据插值求出，计算公式如下所示：其中ugE，ugA，ugD分别为E点的纬向地转风速，A点的纬向地转风速和D点的纬向地转风速，vgE，vgA，vgD分别为E点的经向地转风速，A点的经向地转风速和D 点的经向地转风速；G点的纬向地转风风速和经向地转风速由F、H两点的数据插值求出，计算公式如下所示：其中ugG，ugF，ugH分别为G点的纬向地转风速，F点的纬向地转风速和H点的纬向地转风速，vgG，vgF，vgH分别为G点的经向地转风速，F点的经向地转风速和H点的经向地转风速；B点的纬向地转风风速和经向地转风速由G、E两点的数据插值求出，计算公式如下所示：其中ugB，ugG，ugE分别为B点的纬向地转风速，G点的纬向地转风速和E点的纬向地转风速，vgB ， vgG ， vgE 分别为B点的经向地转风速，G点的经向地转风速和E点的经向地转风速；计算C点纬向地转风风速和经向地转风风速：利用A，I，J，K四个网格点进行插值，得到C点的纬向地转风速ugc和C点的经向地转风速vgc，公式如下所示：其中，ugL，ugA，ugK，ugM，ugJ，ugI，ugC分别为L的纬向地转风风速、A的纬向地转风风速、K的纬向地转风风速、M的纬向地转风风速、J的纬向地转风风速、I的纬向地转风风速和C的纬向地转风风速，vgL，υgA，vgK，vgM，υgJ，vgI，vgC分别为L的经向地转风风速、A的经向地转风风速、K的经向地转风风速、M的经向地转风风速、J的经向地转风风速、I的经向地转风风速和C的经向地转风风速；根据B，C两点的经向和纬向的地转风风速，分别计算B，C两点的地转风风向角度，并根据经向、纬向地转风风速的正负，映射至角度坐标系中，即将纬向地转风风速ug与经向地转风风速vg相除，对计算结果做反正切运算，得到地转风风向角度，再根据经向、纬向的地转风风速的正负，将地转风风向角度映射至角度坐标系中；步骤1‑3‑1‑2，计算网格点A的曲率值：在A点的地转风风速方向上对B，C两点的地转风风向角度进行中央差分运算，其结果即为网格点A的曲率值，该曲率值视为过A点的位势高度等值线在该点处的曲率值，计算公式如下所示：CA＝θC‑θB，其中，CA为A点处的曲率值，θB为B点处的地转风风向角度，θC为C点处的地转风风向角度；步骤1‑3‑1‑3，计算非边界网格点曲率值：遍历位势高度场网格数据中的非边界数据网格点，即行标大于1且小于网格数据行数Row，列标大于1且小于网格数据列数Rank的数据网格点，对其进行如步骤1‑3‑1‑1和步骤1‑3‑1‑2中所述的运算，得到所有非边界数据网格点对应的曲率值；步骤1‑3‑2包括以下步骤：步骤1‑3‑2‑1，插值计算邻近一网格点的地转风风向角度：网格点距为L，其中，O，R，S，T为网格点，O为待求曲率值的网格点，S、R、T分别为O点下方网格点、左下方网格点和左方网格点，OP的长度为β为O点的地转风风向角度，N、Q分别为OT上的点和SR上的点，N点的纬向地转风风速和经向地转风速由O、T两点的数据插值求出，计算公式如下所示：其中，ugN，ugO，ugT分别为N点的纬向地转风速，O点的纬向地转风速和T点的纬向地转风速，vgN，vgO，vgT分别为N点的经向地转风速，O点的经向地转风速和T点的经向地转风速；Q点的纬向地转风风速和经向地转风风速由S、R两点的数据插值求出，计算公式如下所示：其中，ugQ，ugS，ugR分别为Q点的纬向地转风速，S点的纬向地转风速和R点的纬向地转风速，vgQ，vgS，vgR分别为Q点的经向地转风速，S点的经向地转风速和R点的经向地转风速；P点的纬向地转风风速和经向地转风速由N、Q两点的数据插值求出，计算公式如下所示：其中ugP，ugN，ugQ分别为P点的纬向地转风速，N点的纬向地转风速和Q点的纬向地转风速，vgP，vgN，vgQ分别为P点的经向地转风速，N点的经向地转风速和Q点的经向地转风速；根据P点的经向和纬向的地转风风速，计算P点的地转风风向角度，并根据经向、纬向地转风风速的正负，映射至角度坐标系中，即将纬向地转风风速ug与经向地转风风速vg相除，对计算结果做反正切运算，得到地转风风向角度，再根据经向、纬向的地转风风速的正负，将地转风风向角度映射至角度坐标系中；步骤1‑3‑2‑2，计算网格点O的曲率值：在O点的地转风风速方向上对O，P两点的地转风风向角度进行差分运算，其结果即为网格点O的曲率值，该曲率值视为过O点的位势高度等值线在该点处的曲率值，计算公式如下所示：CO＝2*(θP‑θO)，其中，CO为O点处的曲率值，θO为O点处的地转风风向角度，θP为P点处的地转风风向角度；步骤1‑3‑2‑3，计算边界网格点曲率值：遍历网格数据中的边界数据网格点，即行标等于1或Row，或者列标等于1或Rank的数据网格点，对其进行如步骤1‑3‑2‑1和步骤1‑3‑2‑2中所述的运算，得到所有边界数据网格点对应的曲率值；步骤1‑4包括以下步骤：步骤1‑4‑1，比较曲率值大小：遍历位势高度场网格数据中的所有非边界数据网格点，即行标大于1且小于网格数据行数Row，列标大于1且小于网格数据列数Rank的数据网格点，将其与周围八个邻接点的曲率值进行比较，其中，X为非边界数据网格点，a、b、c、d、e、f、g、h为X的顺时针排列的八个邻接点，将X的曲率值与这八个点逐一比较；步骤1‑4‑2，存储曲率值较小邻接点数量：将步骤1‑4‑1中曲率值比X点小的邻接点数量存储起来，用N来表示，遍历所有非边界网格点，得到每一点所对应的的曲率值较小邻接点数量Nij，其中i，j分别为行索引和列索引；步骤2包括以下步骤：步骤2‑1，提取槽线追踪起点：遍历所有非边界网格点的曲率值较小邻接点数量Nij，提取所有Nij＝8的网格点作为槽线追踪起点，存储在集合StartPoint中，并建立已追踪网格点集合MarkedPoint；步骤2‑2，槽线追踪：根据各非边界网格点的曲率值较小邻接点数量Nij，以及槽线追踪角度的约束，对槽线进行两个方向的追踪；其中，步骤2‑2包括以下步骤：步骤2‑2‑1，分类标记网格点，遍历所有非边界网格点，将曲率值较小邻接点数量Nij≤3的网格点的Nij统一设置为0；将Nij≥6的网格点的Nij保持实际值不变；对于Nij＝4或Nij＝5的网格点，如满足如下条件，则该点的Nij仍保持实际值不变，否则即将其Nij设置为0：对于Nij＝5的网格点，假设其在a、b处各有一个曲率值较大的邻接点，则如其第三个较大邻接点在d或e或f或g点，则将该点的Nij仍保持实际值不变，这里的a、b、d、e、f、g均为相对位置，即如该网格点在a、h处各有一个曲率值较大的邻接点，则第三个较大邻接点应在c、d、e、f处，如不满足该条件，则将该点的Nij设置为0；对于Nij＝4的网格点，假设其在a、b处各有一个曲率值比该网格点大的邻接点，如另一对曲率值较大邻接点在d、e点或e、f点或f、g点，则将该点的Nij仍保持实际值不变，这里的a、b、d、e、f、g均为相对位置，即如该网格点在a、h处各有一个曲率值较大的邻接点，则另一对曲率值较大邻接点应在c、d点或d、e点或e、f点，如不满足该条件，则将该点的Nij设置为0；步骤2‑2‑2，选择槽线追踪起点：从满足Nij＝8的点集中按序任意选取一个点作为槽线追踪起点，并建立包含该追踪起点的槽点集合PointListk，其中k是槽点集的索引，k取值为自然数；步骤2‑2‑3，根据槽线追踪起点周围邻接点的Nstartpoint值确定下一追踪槽点，NStartPoint表示槽线追踪起点的曲率值较小邻接点数量；在Nij≠0的网格点的集合中确定下一追踪槽点，具体遵循以下规则：优先选择槽线追踪起点周围Nij最大的邻接点作为下一追踪槽点，如存在两个以上邻接点的Nij值为并列最大，则优先选择其中曲率值更大的邻接点作为下一追踪槽点，并将槽线追踪起点和已追踪槽点加入集合MarkedPoint；步骤2‑2‑4，根据上一追踪槽点，结合约束条件确定下一追踪槽点，直至当前槽线追踪结束；在根据上一追踪槽点确定下一追踪槽点时，应满足以下约束条件：追踪角度偏转最大不超过90°，如上一追踪方向是从a至X点，则下一追踪槽点应从c、d、e、f、g五点中选择，此处的a、X、c、d、e、f、g均为相对位置，在满足追踪角度要求的X的邻接点中，优先选择Nij最大的邻接点作为下一追踪槽点，如存在两个以上邻接点的Nij值为并列最大，则优先选择其中曲率值更大的邻接点作为下一追踪槽点，并将已追踪槽点加入集合MarkedPoint；槽点追踪如满足以下任一约束条件，则当前槽线追踪结束：追踪过程中遇到已标记网格点，即集合MarkedPoint中点，或者由当前点追踪下一槽点时，当前点的八个邻接点Nij值均为0；步骤2‑2‑5，存储已追踪槽线：对已追踪完成的槽线中的槽点存储至集合PointListk中，其中k是槽点集的索引；步骤2‑2‑6，按序遍历槽线追踪起点集合StartPoint，重复步骤2‑2‑3～2‑2‑5直至最后一个槽线追踪起点，最终得到k个原始槽线集合。