CN116027459B - 基于数值天气预报数据的大气加权平均温度的计算方法 - Google Patents
基于数值天气预报数据的大气加权平均温度的计算方法 Download PDFInfo
- Publication number
- CN116027459B CN116027459B CN202211720316.9A CN202211720316A CN116027459B CN 116027459 B CN116027459 B CN 116027459B CN 202211720316 A CN202211720316 A CN 202211720316A CN 116027459 B CN116027459 B CN 116027459B
- Authority
- CN
- China
- Prior art keywords
- grapes
- weighted average
- average temperature
- interpolation
- deviation
- 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
Links
- 238000004364 calculation method Methods 0.000 title abstract description 16
- 238000012417 linear regression Methods 0.000 claims abstract description 26
- 238000000034 method Methods 0.000 claims abstract description 20
- 230000009897 systematic effect Effects 0.000 claims abstract description 16
- 238000012937 correction Methods 0.000 claims abstract description 10
- 241000219094 Vitaceae Species 0.000 claims description 83
- 235000021021 grapes Nutrition 0.000 claims description 83
- 230000002159 abnormal effect Effects 0.000 claims description 17
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000011160 research Methods 0.000 claims description 6
- 230000010354 integration Effects 0.000 claims description 4
- 238000012353 t test Methods 0.000 claims description 4
- 238000000692 Student's t-test Methods 0.000 claims description 2
- 230000008859 change Effects 0.000 abstract description 6
- 230000007547 defect Effects 0.000 abstract description 4
- 230000002123 temporal effect Effects 0.000 abstract 1
- 238000012360 testing method Methods 0.000 description 7
- 230000008901 benefit Effects 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 102100033814 Alanine aminotransferase 2 Human genes 0.000 description 1
- 101000779415 Homo sapiens Alanine aminotransferase 2 Proteins 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 239000005442 atmospheric precipitation Substances 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Complex Calculations (AREA)
Abstract
本发明属于气象温度计算技术领域，具体涉及基于数值天气预报数据的大气加权平均温度的计算方法。本发明基于GRAPES_MESO预报数据和探空站的分层数据分别求取Tm，并且判断两者是否存在系统偏差，基于线性回归模型和球冠谐模型构建系统偏差改正模型，进而输出空间分辨率为0.1°×0.1°、时间分辨率为3 h的Tm格网预报产品，再将格网预报产品结合空间插值和三次样条插值即可实现至少提前24小时预报中国地区的Tm。相比于国际公认经验模型GPT2w具有更高的精度，同时该计算方法还可以较好地感知Tm天气尺度的变化，弥补了其他经验模型的不足。
Description
技术领域
本发明属于气象温度计算技术领域，具体涉及基于数值天气预报数据的大气加权平均温度的计算方法。
背景技术
GNSS信号穿越大气层时由于折射的影响，往往会出现一种被称为对流层延迟的信号弯曲或延迟现象。其中的湿分量可以通过大气加权平均温度(Tm)线性转化为大气可降水量(Precipitation water vapor,PWV)。而大气水汽含量的时空变迁是监测和预报不同尺度天气、气候现象的关键因子。因此，基于GNSS反演大气水汽含量受到国内外学者的广泛关注，进而延伸了出了一个新的学科——GNSS气象学。综上，精确计算Tm是利用GNSS准确反演PWV的重要保障，进而提升基于大气水汽监测和预报不同天气、气候的性能。精确的Tm虽可通过测站的实测气象数据计算获得，但不是所有地区都有实测气象数据。为此，学者们构建了一系列只需测站时空信息驱动的Tm经验模型，这些模型对实时GNSS应用起到了极大的促进作用。然而，多数经验模型是基于Tm多年的平均值、年周期、半年周期和日周期变化振幅构建的，同时假设这些变化是每年重复出现的。故经验模型仅能表征Tm的季节的变化，而无法感知Tm天气尺度的复杂变化，精度有待进一步提升。数值天气预报(Numerical WeatherPrediction,NWP)数据的出现且持续更新，为预报Tm提供了一种新思路。欧洲中尺度天气预报中心(European Centre for Medium-Range Weather Forecasts,ECMWF)和美国国家环境预报中心(National Centers for Environmental Prediction,NCEP)两个常用数值天气预报数据被学者们广泛接受。其中，计算“维也纳映射函数”预测数据——VMF1-FC使用的就是ECMWF预测数据。已有学者研究发现，NWP数据在预测对流层延迟的效果比经验模型更好，因为数值天气预报数据总是基于最近的数据模拟对流层的状态。然而遗憾的是，并不是所有的用户都可以免费获取ECMWF和NCEP数据，尤其对于中国地区的用户。因此，为解决经验模型难以捕捉Tm天气尺度变化的缺点、中国用户难以获得ECMWF和NCEP的难点，本发明利用中国气象局中尺度版本的全球/区域同化和预报系统(GRAPES_MESO)预测数据，开发了一种中国地区预报Tm的计算方法(CTm_FC)，此计算方法能使中国地区用户预测Tm更准确、更便捷。
发明内容
为了解决上述问题，本发明提供了基于数值天气预报数据的大气加权平均温度的计算方法，以解决经验模型无法感知Tm在天气尺度变化的问题，提高中国用户预报Tm的精度和便捷性。具体技术方案如下：
基于数值天气预报数据的大气加权平均温度的计算方法，包括以下步骤：步骤S1，选取研究区域的GRAPES_MESO预报数据计算大气加权平均温度Tm，记为GRAPES Tm；
步骤S2，选取与步骤S1的GRAPES_MESO预报数据处于相同时间段的M个探空站的分层数据积分求取大气加权平均温度Tm，记为RS Tm；
步骤S3，对步骤S1中计算得到的GRAPES Tm进行水平插值和垂直插值得到任意位置与任意高度的大气加权平均温度Tm；
步骤S4，将步骤S3中对应M个探空站所在位置插值后的GRAPES Tm与RS Tm进行比较，剔除异常值后计算GRAPES Tm与RS Tm的差值的平均偏差、标准偏差、均方根误差，判断是否需要改正GRAPES Tm与RS Tm间的系统偏差，如若需要，则进入步骤S5；
步骤S5，选择线性回归模型，将RS Tm作为参考值，针对对应探空站的GRAPES Tm进行系统偏差改正，得到各探空站的线性回归模型的比例偏差参数a和恒定偏差参数b；
步骤S6，构建球冠谐模型，以步骤S5得到的各探空站的线性回归模型的比例偏差参数a和恒定偏差参数b确定球冠谐模型的阶数、次数以及相关系数，进而构建合适大气加权平均温度Tm的比例偏差参数a和恒定偏差参数b展开的球冠谐模型；
步骤S7，使用已经确定的球冠谐模型计算各个GRAPES网格点的线性回归模型的比例偏差参数a、恒定偏差参数b，并修正GRAPES Tm，进而获得大气加权平均温度Tm的格网预报产品，记为CTm_FC；
步骤S8，对格网预报产品值CTm_FC进行空间插值和三次样条插值，得到至少提前24小时预报中国地区的大气加权平均温度Tm值，记为CTm_FC。
优选地，所述GRAPES_MESO预报数据包括提前24小时的气压、温度、水汽压和位势高数据。
优选地，所述步骤S1和步骤S2采用以下公式计算大气加权平均温度Tm：
其中，ei和Ti分别表示i层水汽压和温度，e指水汽压、T代表温度、h为高程，Δhi指i层的高度。
优选地，所述步骤S3中的水平插值是是将离目标位置最近的四个格网节点的GRAPES Tm向目标位置进行插值，得到目标位置的GRAPES Tm，记为Tmt；具体包括以下步骤：
设离目标位置最近的四个网格节点的GRAPES Tm分别为：Tm1、Tm2、Tm3、Tm4，其中GRAPES Tm1、GRAPES Tm2、GRAPES Tm3、GRAPES Tm4按逆时针顺序排列；计算待插值点在水平面内的横向和纵向的权重，具体如下：
其中S1为在插值水平面内格网点Tm1或Tm4在纵向上到待插值点Tmt的距离；S2为在插值水平面内格网点Tm2或Tm3在纵向上到待插值点Tmt的距离；S3为在插值水平面内格网点Tm1或Tm2在横向上到待插值点Tmt的距离；S4为在插值水平面内格网点Tm4或Tm3在横向上到待插值点Tmt的距离；P1为S1在插值水平面内纵向上的权重；P2为S2在插值水平面内纵向上的权重；P3为S3在插值水平面内横向上的权重；P4为S4在插值水平面内横向上的权重；
计算目标位置的待插值的Tmt，具体如下：
Tmt＝P3×(Tm1×P1+Tm2×P2)+P4×(Tm4×P1+Tm3×P2)； (6)。
优选地，所述步骤S3采用二阶傅立叶函数进行垂直插值，具体如下：
Tmt＝Tm0+a1·(cos(w×h0)-cos(w×ht))+b1·(sin(w×h0)-sin(w×ht))+a2·(cos(2×w×h0)-cos(2×w×ht))+b2·(sin(2×w×h0)-sin(2×w×ht))； (7)
其中，Tmt表示目标高度ht的Tm，Tm0表示格网点高度h0处的Tm，a1、b1、a2、b2和w是通过公式(7)对Tm垂直廓线进行拟合求解所得二阶傅立叶系数。
优选地，所述步骤S4中具体包括以下步骤：
步骤S41，去除M个样本对中的负值；
步骤S42，计算GRAPES Tm与RS Tm的差值；
步骤S43，检验所述差值是否服从正态分布并且计算差值的标准偏差；
步骤S44，以标准偏差的3倍作为异常阈值，若对应GRAPES Tm与RS Tm的差值大于异常阈值，则判定对应样本异常并进行剔除；
步骤S45，剔除异常值后计算GRAPES Tm与RS Tm的Tm差值的平均偏差、标准偏差、均方根误差；
步骤S46，判断Tm差值的平均偏差的绝对值是否大于均方根误差的40％-50％，若是，则判断需要改正GRAPES Tm与RS Tm间的系统偏差。
优选地，所述步骤S5中具体包括以下步骤：
步骤S51，应用线性回归模型在M个探空站对GRAPES Tm和RS Tm利用公式(8)拟合和求解每个探空站点的线性回归模型的比例偏差参数a和恒定偏差参数b；
Yi＝ai×Xi+bi； (8)
其中，i为不同探空站点的序号，i＝1,2，…，M，Xi代表第i个探空站的GRAPES Tm矩阵，Yi代表对应的第i个探空站的RS Tm矩阵，ai为第i个探空站的线性回归模型的比例偏差参数，bi为第i个探空站的线性回归模型的恒定偏差参数；
步骤S52，使用t检验检查步骤S1求得的每个探空站点的参数a、b的显著性；
步骤S53，对参数a的显著性检查中超过90％的站点p值小于0.05，以及参数b的显著性检查中超过50％的站点p值小于0.05，则输出最终的各个探空站点的参数a、b。
优选地，所述步骤S6中具体包括以下步骤：
步骤S61，构建球冠谐模型，所述球冠谐模型具体如下：
其中，C代表比例偏差参数a或恒定偏差参数b，ae是地球半径，r、θ、λ分别代表地表某点在球面坐标系下的向径、余纬和经度，是缔合勒让德函数，nk(m)和m分别表示球冠谐模型的阶数和次数，下标k是用于标识nk(m)顺序的整数，m是整数的实数；nk(m)是实数；N表示a或b的最大展开阶次；/>和/>是球冠谐模型的待定系数；其中阶数nk(m)和次数m相等；
步骤S62，在研究区域外设置虚拟参考站并将虚拟参考站中的参数a、b分别设置为1和0；
步骤S63，采用以阶数递增的球冠谐模型来分别对参数a和b进行拟合，并以拟合得到的均方根误差来评估拟合精度，进而得到球冠谐模型合适的阶数与次数；步骤S64，采用GRAPES Tm的参数a和b值来解算球冠谐模型的待定系数，进而构建合适大气加权平均温度Tm的参数a和b展开的球冠谐模型。
优选地，用于构建球冠谐模型的参数a和b需满足以下条件：
(1)用于计算参数a、b的GRAPES Tm与RS Tm数据集间相关系数大于等于0.8；(2)单站求解所得a、b超出对应所有参数的三倍标准偏差，视为异常值剔除。
本发明的有益效果为：本发明提供了基于数值天气预报数据的大气加权平均温度的计算方法，相比于国际公认经验模型GPT2w更具有更高的精度，同时该计算方法还可以较好地感知Tm天气尺度的变化，弥补了其他经验模型的不足。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍。在所有附图中，类似的元件或部分一般由类似的附图标记标识。附图中，各元件或部分并不一定按照实际的比例绘制。
图1为本发明的流程示意图；
图2为本发明空间插值示意图；
图3为本发明拟合所得Tm的系统差改正系数a和b的显著性检验p值图；
图4为本发明BAGUIO和LUCKNOW站改正前后GRAPES Tm和RS Tm时间序列对比图；
图5为本发明GRAPES Tm的a和b拟合RMS与SCH模型阶数的对应关系图；
图6为本发明CTm_FC相对GPT2w的RMS降低率分布图；
图7为本发明CTm_FC和GPT2w模型预报Tm的bias和RMS功率谱图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
应当理解，当在本说明书和所附权利要求书中使用时，术语“包括”和“包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
还应当理解，在本发明说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本发明。如在本发明说明书和所附权利要求书中所使用的那样，除非上下文清楚地指明其它情况，否则单数形式的“一”、“一个”及“该”意在包括复数形式。
还应当进一步理解，在本发明说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
如图1所示，本发明的具体实施方式提供了基于数值天气预报数据的大气加权平均温度的计算方法，包括以下步骤：
步骤S1，选取研究区域的GRAPES_MESO预报数据计算大气加权平均温度Tm，记为GRAPES Tm。所述GRAPES_MESO预报数据包括提前预报24小时的气压、温度、水汽压和位势高数据。
步骤S2，选取与步骤S1的GRAPES_MESO预报数据处于相同时间段的M个探空站的分层数据积分求取大气加权平均温度Tm，记为RS Tm。
在本实施例中，以中国及周边地区为例，选择研究区域(15°N-55°N和70°E-135°E)。本发明选取2017年GRAPES_MESO往前预报24小时的气压、温度、水汽压和位势高数据计算格网Tm预报产品，数据从中国气象数据网下载。同时选取2017年146个探空站的分层数据积分求取Tm，用作检验。探空站的分层数据的数据来源为Integrated Global RadiosondeArchive(IGRA)，其时间分辨率为12h。
GRAPES_MESO预报数据与探空站的分层数据均提供了垂直分层的气压、温度和相对湿度等气象参数数据。因此，本发明将采用以下公式计算大气加权平均温度Tm：
其中，ei和Ti分别表示i层水汽压和温度，e指水汽压、T代表温度、h为高程，Δhi指i层的高度。
步骤S3，如图2所示，对步骤S1中计算得到的GRAPES Tm进行水平插值和垂直插值得到任意位置与任意高度的大气加权平均温度Tm。
参考图2(b)，本发明中水平插值采用双线性插值法，水平插值是将离目标位置最近的四个格网节点的GRAPES Tm向目标位置进行插值，得到目标位置的GRAPES Tm，记为Tmt；具体包括以下步骤：
设离目标位置最近的四个网格节点的GRAPES Tm分别为：Tm1、Tm2、Tm3、Tm4，其中GRAPES Tm1、GRAPES Tm2、GRAPES Tm3、GRAPES Tm4按逆时针顺序排列；计算待插值点在水平面内的横向和纵向的权重，具体如下：
其中S1为在插值水平面内格网点Tm1或Tm4在纵向上到待插值点Tmt的距离；S2为在插值水平面内格网点Tm2或Tm3在纵向上到待插值点Tmt的距离；S3为在插值水平面内格网点Tm1或Tm2在横向上到待插值点Tmt的距离；S4为在插值水平面内格网点Tm4或Tm3在横向上到待插值点Tmt的距离；P1为S1在插值水平面内纵向上的权重；P2为S2在插值水平面内纵向上的权重；P3为S3在插值水平面内横向上的权重；P4为S4在插值水平面内横向上的权重；
计算目标位置的待插值的Tmt，具体如下：
Tmt＝P3×(Tm1×P1+Tm2×P2)+P4×(Tm4×P1+Tm3×P2)； (6)。
如图2(a)所示，对于垂直插值，具体是采用二阶傅立叶函数进行垂直插值，具体如下：
Tmt＝Tm0+a1·(cos(w×h0)-cos(w×ht))+b1·(sin(w×h0)-sin(w×ht))+a2·(cos(2×w×h0)-cos(2×w×ht))+b2·(sin(2×w×h0)-sin(2×w×ht))； (7)
其中，Tmt表示目标高度ht的Tm，Tm0表示格网点高度h0处的Tm，a1、b1、a2、b2和w是通过公式(7)对Tm垂直廓线进行拟合求解所得二阶傅立叶系数。
步骤S4，将步骤S3中对应M个探空站所在的位置的插值后的GRAPES Tm与RS Tm进行比较，剔除异常值后计算GRAPES Tm与RS Tm的Tm差值的平均偏差、标准偏差、均方根误差，判断是否需要改正GRAPES Tm与RS Tm间的系统偏差，如若需要，则进入步骤S5。
因探空数据具有可靠的准确性，所以被大量应用于检测其他产品的可靠性。因此，本发明以RS Tm为参考验证插值后的GRAPES Tm。本发明对往前预报24小时的00:00UTC和12:00UTC的GRAPES Tm插值到146个探空站中，并与RS Tm做比较。具体包括以下步骤：
步骤S41，去除146个样本对中的负值；
步骤S42，计算GRAPES Tm与RS Tm的差值；
步骤S43，检验所述差值是否服从正态分布并且计算差值的标准偏差；
步骤S44，以标准偏差的3倍作为异常阈值，若对应GRAPES Tm与RS Tm的差值大于异常阈值，则判定对应样本异常并进行剔除；
步骤S45，剔除异常值后计算GRAPES Tm与RS Tm的Tm差值的平均偏差(Bias)、标准偏差(Standard Deviation,STD)、均方根误差(Root mean square,RMS)。结果见表1。
表1探空数据验证GRAPES预测Tm的精度
步骤S46，判断Tm差值的平均偏差的绝对值是否大于均方根误差的40％-50％，若是，则判断需要改正GRAPES Tm与RS Tm间的系统偏差。
由表1可知，GRAPES Tm的bias和RMS为-1.0K和1.9K，由于bias的绝对值已经超过RMS的50％。这种量级的bias对RMS影响较大，可视为系统偏差，需要对其进行合理修正。
步骤S5，选择线性回归模型，将RS Tm作为参考值，针对对应探空站的GRAPES Tm进行系统偏差改正，得到各探空站的线性回归模型的比例偏差参数a和恒定偏差参数b。具体包括以下步骤：
步骤S51，应用线性回归模型在146个探空站对GRAPES Tm和RS Tm利用公式(8)拟合和求解每个探空站点的线性回归模型的比例偏差参数a和恒定偏差参数b；
Yi＝ai×Xi+bi； (8)
其中，i为不同探空站点的序号，i＝1,2，...，M，Xi代表第i个探空站的GRAPES Tm矩阵，Yi代表对应的第i个探空站的RS Tm矩阵，ai为第i个探空站的线性回归模型的比例偏差参数，bi为第i个探空站的线性回归模型的恒定偏差参数；
步骤S52，使用t检验检查步骤S1求得的每个探空站点的参数a、b的显著性，t检验基于95％置信度展开；
步骤S53，对参数a的显著性检查中超过90％的站点p值小于0.05，以及参数b的显著性检查中超过50％的站点p值小于0.05，则输出最终的各个探空站点的参数a、b。
由图3结果显示，对参数a的显著性检查中超过90％的站点p值小于0.05，b中超过50％的站点p值小于0.05，说明公式(8)和解算所得a、b具有统计意义。利用解算所得a、b将GRAPES Tm改正后，bias有显著改进并接近零。图4给出了经改正后的两个示例站的GRAPESTm与RS Tm时间序列，显然经系统偏差改正后的GRAPES Tm与RS Tm吻合更好，印证了本发明使用线性回归模型对GRAPES Tm的系统偏差改正具有良好的可靠性与准确性的。
使用线性回归模型仅可获得探空站单站的GRAPES Tm系统差改正，为获取研究区内任何位置的GRAPES Tm系统偏差，本发明使用球冠谐模型(SCH)来拟合参数a、b与空间的关系，具体如以下步骤。
步骤S6，构建球冠谐模型，以步骤S5得到的各探空站的线性回归模型的比例偏差参数a和恒定偏差参数b确定球冠谐模型的阶数、次数以及相关系数，进而构建合适大气加权平均温度Tm的比例偏差参数a和恒定偏差参数b展开的球冠谐模型。
具体包括以下步骤：
步骤S61，构建球冠谐模型，所述球冠谐模型具体如下：
其中，C代表比例偏差参数a或恒定偏差参数b，ae是地球半径，r、θ、λ分别代表地表某点在球面坐标系下的向径、余纬和经度，是缔合勒让德函数，nk(m)和m分别表示球冠谐模型的阶数和次数，下标k是用于标识nk(m)顺序的整数，m是整数的实数；nk(m)是实数；N表示a或b的最大展开阶次；/>和/>是球冠谐模型的待定系数；其中阶数nk(m)和次数m相等。
步骤S62，为了防止边缘效应，在研究区域外设置虚拟参考站并将虚拟参考站中的参数a、b分别设置为1和0。为了确保构建出的模型具有一定准确性，关键在于解算模型采用的高质量的a、b数据。因此用于构建球冠谐模型的参数a和b需满足以下条件：
(1)用于计算参数a、b的GRAPES Tm与RS Tm数据集间相关系数大于等于0.8；(2)单站求解所得a、b超出对应所有参数的三倍标准偏差，视为异常值剔除。步骤S63，为了获得球冠谐模型的最佳阶数与次数，采用以阶数递增的球冠谐模型来分别对参数a和b进行拟合，并以拟合得到的均方根误差来评估拟合精度，进而得到球冠谐模型合适的阶数与次数。
由图5知，GRAPES Tm的参数a和b阶数最佳可能为10或11，此时均方根误差最小。但将参数a和b的阶数确定为10进行拟合时，发现西北地区以及青藏高原地区出现多处畸变区域；继续降低阶数到9时重新拟合，青藏高原地区仍然存在明显畸变；阶数降低为8时，在西北地区依然存在畸变区域；直至阶数降低到7时进行拟合，整个研究区域内已经不存在明显畸变区域。出现畸变区域时代表该阶数下参数a、b出现了过度拟合的现象，会对最终结果造成影响。所以为保证整体精度较高且过拟合现象尽可能少，最终本发明将SCH模型的阶数定为7来拟合GRAPES Tm的参数a、b。
步骤S64，阶数和次数确定后，采用GRAPES Tm的参数a和b值来解算球冠谐模型的待定系数，进而构建合适大气加权平均温度Tm的参数a和b展开的球冠谐模型。
步骤S7，使用已经确定的球冠谐模型计算各个GRAPES网格点的线性回归模型的比例偏差参数a、恒定偏差参数b，并修正GRAPES Tm，进而获得大气加权平均温度Tm的格网预报产品，记为CTm_FC。
步骤S8，对格网预报产品值CTm_FC进行空间插值和三次样条插值，得到至少提前24小时预报中国地区的大气加权平均温度Tm值，记为CTm_FC。
本发明利用GRAPES_MESO预测产品计算Tm，发现其GRAPES Tm与RS Tm存在不可忽略的系统偏差，基于线性回归模型和球冠谐模型构建系统偏差改正模型，进而输出空间分辨率为0.1°×0.1°、时间分辨率为3h的Tm格网预报产品，再将格网预报产品结合空间插值和三次样条插值即可实现至少提前24小时预报中国地区的Tm。
本发明相比于国际公认经验模型GPT2w更具有更高的精度，同时该计算方法还可以较好地感知Tm天气尺度的变化，弥补了其他经验模型的不足。本发明基于三个指标评价本发明的有效性，指标为偏差(Bias)、标准偏差(Standard Deviation,STD)、均方根误差(Root mean square,RMS)。表2给出了CTm_FC和GPT2w的bias、STD和RMS，显然，CTm_FCbias、STD与RMS比GPT2w均有明显优化，CTm_FC的bias降低了85.8％；用RMS来衡量，CTm_FC总体精度提高了60.5％。上述结果证明，CTm_FC预报Tm比GPT2w具有更高的精度，这种优势来源于CTm_FC可以敏锐地感知Tm天气尺度的变化。
表2探空数据检验CTm_FC和GPT2w的bias、STD和RMS
本发明还选取了2017年146个探空站的RS Tm对计算方法的结果进行验证，图6是CTm_FC相对GPT2w在各站的RMS提升率。可以明显看到，除极个别站外CTm_FC的RMS均有明显的减小，RMS平均降低率为57.7％，精度明显高于GPT2w。
为了探讨本发明提出的CTm_FC相对于GPT2w精度有所提高的原因，本发明采用FFT对日均bias和RMS进行分析。图7为GPT2w和CTm_FC日均bias和RMS的功率谱。显然，根据CTm_FC的功率谱可以观察到，GPT2w的bias和RMS在1-10天的周期上都具有很大的幅值，这是因为其未模型化Tm天气尺度的变化造成的，然而CTm_FC bias和RMS在相同的时间尺度上，幅值显著减小，这表明CTm_FC很好地捕捉到了Tm的天气尺度变化。这是CTm_FC相对GPT2w等模型的另一突出优势，也是精度优势重要来源。
本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。
在本申请所提供的实施例中，应该理解到，单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元可结合为一个单元，一个单元可拆分为多个单元，或一些特征可以忽略等。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。
Claims (8)
1.基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，包括以下步骤：
步骤S1，选取研究区域的GRAPES_MESO预报数据计算大气加权平均温度Tm，记为GRAPESTm；
步骤S2，选取与步骤S1的GRAPES_MESO预报数据处于相同时间段的M个探空站的分层数据积分求取大气加权平均温度Tm，记为RS Tm；
步骤S3，对步骤S1中计算得到的GRAPES Tm进行水平插值和垂直插值得到任意位置的大气加权平均温度Tm；
步骤S4，将步骤S3中对应M个探空站所在的位置的插值后的GRAPES Tm与RS Tm进行比较，剔除异常值后计算GRAPES Tm与RS Tm的Tm差值的平均偏差、标准偏差、均方根误差，判断是否需要改正GRAPES Tm与RS Tm间的系统偏差，如若需要，则进入步骤S5；具体包括以下步骤：
步骤S41，去除M个样本对中的负值；
步骤S42，计算GRAPES Tm与RS Tm的差值；
步骤S43，检验所述差值是否服从正态分布并且计算差值的标准偏差；
步骤S44，以标准偏差的3倍作为异常阈值，若对应GRAPES Tm与RS Tm的差值大于异常阈值，则判定对应样本异常并进行剔除；
步骤S45，剔除异常值后计算GRAPES Tm与RS Tm的Tm差值的平均偏差、标准偏差、均方根误差；
步骤S46，判断Tm差值的平均偏差的绝对值是否大于均方根误差的40％-50％，若是，则判断需要改正GRAPES Tm与RS Tm间的系统偏差；
步骤S5，选择线性回归模型，将RS Tm作为参考值，针对对应探空站的GRAPES Tm进行系统偏差改正，得到各探空站的线性回归模型的比例偏差参数a和恒定偏差参数b；
步骤S6，构建球冠谐模型，以步骤S5得到的各探空站的线性回归模型的比例偏差参数a和恒定偏差参数b确定球冠谐模型的阶数、次数以及相关系数，进而构建合适大气加权平均温度Tm的比例偏差参数a和恒定偏差参数b展开的球冠谐模型；
步骤S7，使用已经确定的球冠谐模型计算各个GRAPES网格点的线性回归模型的比例偏差参数a、恒定偏差参数b，并修正GRAPES Tm，进而获得大气加权平均温度Tm的格网预报产品，记为CTm_FC；
步骤S8，对格网预报产品CTm_FC进行空间插值和三次样条插值，得到至少提前24小时预报的大气加权平均温度Tm值，记为CTm_FC。
2.根据权利要求1所述的基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，所述GRAPES_MESO预报数据包括提前24小时的气压、温度、水汽压和位势高数据。
3.根据权利要求2所述的基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，所述步骤S1和步骤S2采用以下公式计算大气加权平均温度Tm：
其中，ei和Ti分别表示i层水汽压和温度，e指水汽压、T代表温度、h为高程，Δhi指i层的高度。
4.根据权利要求1所述的基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，所述步骤S3中的水平插值是是将离目标位置最近的四个格网节点的GRAPES Tm向目标位置进行插值，得到目标位置的GRAPES Tm，记为Tmt；具体包括以下步骤：
设离目标位置最近的四个网格节点的GRAPES Tm分别为：Tm1、Tm2、Tm3、Tm4，其中GRAPESTm1、GRAPES Tm2、GRAPES Tm3、GRAPES Tm4按逆时针顺序排列；
计算待插值点在水平面内的横向和纵向的权重，具体如下：
其中S1为在插值水平面内格网点Tm1或Tm4在纵向上到待插值点Tmt的距离；S2为在插值水平面内格网点Tm2或Tm3在纵向上到待插值点Tmt的距离；S3为在插值水平面内格网点Tm1或Tm2在横向上到待插值点Tmt的距离；S4为在插值水平面内格网点Tm4或Tm3在横向上到待插值点Tmt的距离；P1为S1在插值水平面内纵向上的权重；P2为S2在插值水平面内纵向上的权重；P3为S3在插值水平面内横向上的权重；P4为S4在插值水平面内横向上的权重；
计算目标位置的待插值的Tmt，具体如下：
Tmt＝P3×(Tm1×P1+Tm2×P2)+P4×(Tm4×P1+Tm3×P2)； (6)。
5.根据权利要求1所述的基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，所述步骤S3采用二阶傅立叶函数进行垂直插值，具体如下：
其中，Tmt表示目标高度ht的Tm，Tm0表示格网点高度h0处的Tm，a1、b1、a2、b2和w是通过公式(7)对Tm垂直廓线进行拟合求解所得二阶傅立叶系数。
6.根据权利要求1所述的基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，所述步骤S5中具体包括以下步骤：
步骤S51，应用线性回归模型在M个探空站对GRAPES Tm和RS Tm利用公式(8)拟合和求解每个探空站点的线性回归模型的比例偏差参数a和恒定偏差参数b；
Yi＝ai×Xi+bi； (8)
其中，i为不同探空站点的序号，i＝1,2，…，M，Xi代表第i个探空站的GRAPES Tm矩阵，Yi代表对应的第i个探空站的RS Tm矩阵，ai为第i个探空站的线性回归模型的比例偏差参数，bi为第i个探空站的线性回归模型的恒定偏差参数；
步骤S52，使用t检验检查步骤S1求得的每个探空站点的参数a、b的显著性；
步骤S53，对参数a的显著性检查中超过90％的站点参数a的p值小于0.05，以及参数b的显著性检查中超过50％的站点参数b的p值小于0.05，则输出最终的各个探空站点的参数a、b。
7.根据权利要求1所述的基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，所述步骤S6中具体包括以下步骤：
步骤S61，构建球冠谐模型，所述球冠谐模型具体如下：
其中，C代表比例偏差参数a或恒定偏差参数b，ae是地球半径，r、θ、λ分别代表地表某点在球面坐标系下的向径、余纬和经度，是缔合勒让德函数，nk(m)和m分别表示球冠谐模型的阶数和次数，下标k是用于标识nk(m)顺序的整数，m是整数的实数；nk(m)是实数；N表示a或b的最大展开阶次；/>和/>是球冠谐模型的待定系数；其中阶数nk(m)和次数m相等；
步骤S62，在研究区域外设置虚拟参考站并将虚拟参考站中的参数a、b分别设置为1和0；
步骤S63，采用以阶数递增的球冠谐模型来分别对参数a和b进行拟合，并以拟合得到的均方根误差来评估拟合精度，进而得到球冠谐模型合适的阶数与次数；
步骤S64，采用GRAPES Tm的参数a和b值来解算球冠谐模型的待定系数，进而构建合适大气加权平均温度Tm的参数a和b展开的球冠谐模型。
8.根据权利要求7所述的基于数值天气预报数据的大气加权平均温度的计算方法，其特征在于，用于构建球冠谐模型的参数a和b需满足以下条件：
(1)用于计算参数a、b的GRAPES Tm与RS Tm数据集间相关系数大于等于0.8；
(2)单站求解所得a、b超出对应所有参数的三倍标准偏差，视为异常值剔除。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211720316.9A CN116027459B (zh) | 2022-12-30 | 2022-12-30 | 基于数值天气预报数据的大气加权平均温度的计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211720316.9A CN116027459B (zh) | 2022-12-30 | 2022-12-30 | 基于数值天气预报数据的大气加权平均温度的计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116027459A CN116027459A (zh) | 2023-04-28 |
| CN116027459B true CN116027459B (zh) | 2024-05-14 |
Family
ID=86071783
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211720316.9A Active CN116027459B (zh) | 2022-12-30 | 2022-12-30 | 基于数值天气预报数据的大气加权平均温度的计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116027459B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117390875B (zh) * | 2023-10-27 | 2024-03-12 | 长安大学 | 一种大气加权平均温度模型的构建方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3343249A1 (de) * | 2017-01-03 | 2018-07-04 | UBIMET GmbH | Verfahren zur erhöhung der räumlichen auflösung einer wettervorhersage |
| CN108416031A (zh) * | 2018-03-12 | 2018-08-17 | 南京恩瑞特实业有限公司 | Nriet气象多源探测资料融合分析系统 |
| CN110378540A (zh) * | 2019-08-02 | 2019-10-25 | 桂林理工大学 | 一种适用于广西北部湾地区的大气加权平均温度计算方法 |
| WO2021164480A1 (zh) * | 2020-02-17 | 2021-08-26 | 东南大学 | 基于空间位置的加权平均温度计算方法 |
| CN113639893A (zh) * | 2021-06-29 | 2021-11-12 | 东南大学 | 一种基于多气象因子的近地加权平均温度信息获取方法 |
| CN113804318A (zh) * | 2021-10-11 | 2021-12-17 | 南京信息工程大学 | 一种获取加权平均温度的数据融合方法及计算设备 |
| WO2022048694A1 (zh) * | 2021-03-17 | 2022-03-10 | 山东科技大学 | 一种基于球谐展开的gnss单点定位方法 |
-
2022
- 2022-12-30 CN CN202211720316.9A patent/CN116027459B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3343249A1 (de) * | 2017-01-03 | 2018-07-04 | UBIMET GmbH | Verfahren zur erhöhung der räumlichen auflösung einer wettervorhersage |
| CN108416031A (zh) * | 2018-03-12 | 2018-08-17 | 南京恩瑞特实业有限公司 | Nriet气象多源探测资料融合分析系统 |
| CN110378540A (zh) * | 2019-08-02 | 2019-10-25 | 桂林理工大学 | 一种适用于广西北部湾地区的大气加权平均温度计算方法 |
| WO2021164480A1 (zh) * | 2020-02-17 | 2021-08-26 | 东南大学 | 基于空间位置的加权平均温度计算方法 |
| WO2022048694A1 (zh) * | 2021-03-17 | 2022-03-10 | 山东科技大学 | 一种基于球谐展开的gnss单点定位方法 |
| CN113639893A (zh) * | 2021-06-29 | 2021-11-12 | 东南大学 | 一种基于多气象因子的近地加权平均温度信息获取方法 |
| CN113804318A (zh) * | 2021-10-11 | 2021-12-17 | 南京信息工程大学 | 一种获取加权平均温度的数据融合方法及计算设备 |
Non-Patent Citations (1)
| Title |
|---|
| 顾及多因子影响的中国区域Tm模型精化研究;李浩杰, 刘立龙, 黄良珂;《大地测量与地球动力学》;20220430;第42卷(第4期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116027459A (zh) | 2023-04-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104021424B (zh) | 用于预测风场中的风机的输出功率的方法和装置 | |
| US20070244645A1 (en) | Gas-condition predicting device and diffusion-condition predicting system | |
| CN111027175A (zh) | 基于耦合模型集成模拟的洪水对社会经济影响的评估方法 | |
| CN109784563B (zh) | 一种基于虚拟测风塔技术的超短期功率预测方法 | |
| CN116027459B (zh) | 基于数值天气预报数据的大气加权平均温度的计算方法 | |
| CN112035448A (zh) | 一种综合地基gnss水汽和气象要素的神经网络短临降水预报方法 | |
| CN103413036B (zh) | 一种连续化的森林火险天气等级预报方法 | |
| CN110378540A (zh) | 一种适用于广西北部湾地区的大气加权平均温度计算方法 | |
| CN109543907B (zh) | 一种复杂地形风资源评估方法及其装置 | |
| CN114463947B (zh) | 一种基于时空网络卷积模型的对流性致灾强风预警预报方法 | |
| CN116910041B (zh) | 一种基于尺度分析的遥感降水产品的逐日订正方法 | |
| CN116449331B (zh) | 一种基于w波段雷达和气象卫星的沙尘粒子数浓度估算方法 | |
| Guan et al. | Development of verification methodology for extreme weather forecasts | |
| Liu et al. | Research on data correction method of micro air quality detector based on combination of partial least squares and random forest regression | |
| CN115049013A (zh) | 一种联合线性和svm的短时降雨预警模型融合方法 | |
| Sá et al. | Two step calibration method for ozone low-cost sensor: Field experiences with the UrbanSense DCUs | |
| CN109948175B (zh) | 基于气象数据的卫星遥感反照率缺失值反演方法 | |
| CN109583095A (zh) | 基于混合统计动力模型的西北太平洋台风延伸期预报方法 | |
| CN113311509B (zh) | 一种基于神经网络的mwhts对海面气压的敏感性测试方法 | |
| CN113176420B (zh) | 一种针对电网杆塔点的风速预报订正系统 | |
| CN112034007A (zh) | 一种微波辐射计间接测量露点温度方法及系统 | |
| CN112380778A (zh) | 一种基于海温的气象干旱预报方法 | |
| CN112528566A (zh) | 基于AdaBoost训练模型的空气质量数据实时校准方法及系统 | |
| KR20210055975A (ko) | 대기모델 바람예측 자료의 검보정 고도화를 통한 모델 예측정확도 개선을 위한 장치 및 방법 | |
| Sperati et al. | A new Wind Atlas to support the expansion of the Italian wind power fleet |
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 |