2023年全国研究生数学建模竞赛华为杯F题强对流降水临近预报求解全过程文档及程序
2023年全国研究生数学建模竞赛华为杯
F题 强对流降水临近预报
原题再现:
我国地域辽阔,自然条件复杂,因此灾害性天气种类繁多,地区差异大。其中,雷雨大风、冰雹、龙卷、短时强降水等强对流天气是造成经济损失、危害生命安全最严重的一类灾害性天气[1]。以2022年为例,我国强对流天气引发风雹灾害造成的死亡失踪人数和直接经济损失分别占73%和69%。由于强对流天气具有突发性和局地性强、生命史短、灾害重等特点,其短时(0~12小时)和临近(0至2小时)预报通常也是天气预报业务中的难点。
传统强对流天气临近预报主要依靠雷达等观测资料,结合风暴识别、追踪技术进行雷达外推预报,即通过外推的方法得到未来时刻的雷达反射率因子,并进一步使用雷达反射率因子和降水之间的经验性关系(即Z-R关系)估计未来时刻的降水量[2]。近年来,随着大数据的积累和计算机算力的发展,人工智能及深度学习技术发展迅速。深度学习方法是一类数据驱动的方法,理论上其性能随着训练数据量增大而提升,因此很适合有大量雷达观测数据积累的短临预报领域。目前国际上主要有两类基于深度学习的短临预报模型,一类基于卷积神经网络(Convolutional Neural Networks, CNNs),如U-Net等模型[3];另一类基于循环神经网络(Recurrent Neural Networks, RNNs),如ConvLSTM、DGMR等模型[4, 5]。
雨滴在降落过程中受到空气阻力作用,形状可呈扁球形或馒头形,并且一般来说越大的雨滴越扁。因此,雨滴对水平偏振(电场振动方向在水平面内)的电磁波和垂直偏振(电场振动方向在垂直平面内)的电磁波的反射特征是不一样的。传统雷达仅能发射和接收一个偏振方向上的电磁波,而新型的双偏振雷达可同时发射和接收在水平和垂直两个偏振方向的电磁波,可以根据两个偏振方向上的回波的强度差别、相位关系等信息获得降水粒子的大小、相态、含水量等信息[6],这些信息被统称为微物理信息。近年来研究表明,双偏振雷达变量反映的微物理信息里包含了对流系统的演变状态、空间动力结构等关键信息[7, 8]。因此,双偏振雷达变量的应用,理论上对于强对流预报有重要意义。
为了更好地应用双偏振雷达改进强对流降水短临预报,请回答以下问题:
1.如何有效应用双偏振变量改进强对流预报,仍是目前气象预报的重点难点问题。请利用题目提供的数据,建立可提取用于强对流临近预报双偏振雷达资料中微物理特征信息的数学模型。临近预报的输入为前面一小时(10帧)的雷达观测量(ZH 、ZDR、KDP),输出为后续一小时(10帧)的ZH预报。
2.当前一些数据驱动的算法在进行强对流预报时,倾向于生成接近于平均值的预报,即存在“回归到平均(Regression to the mean)”问题,因此预报总是趋于模糊。在问题1的基础上,请设计数学模型以缓解预报的模糊效应,使预报出的雷达回波细节更充分、更真实。
3.请利用题目提供的ZH、ZDR和降水量数据,设计适当的数学模型,利用ZH及ZDR进行定量降水估计。模型输入为ZH和ZDR,输出为降水量。(注意:算法不可使用KDP变量。)
4.请设计数学模型来评估双偏振雷达资料在强对流降水临近预报中的贡献,并优化数据融合策略,以便更好地应对突发性和局地性强的强对流天气。
名词解释:
1.双偏振雷达: 一种新型的气象探测雷达,能够提供比传统雷达更丰富的物理信息。它通过测量降水粒子对水平和垂直两个方向上的电磁波的反射情况,来获取降水粒子的大小、相态、含水量等信息。这些信息被统称为微物理信息,能够帮助我们更好地预测强对流天气。双偏振雷达最常用的三个变量为:1)ZH,水平反射率因子,即水平方向的回波强度,单位通常为dBZ,主要反映降水的强弱;2)ZDR,差分反射率,即水平和垂直方向回波强度的差异,主要反映了观测区域的降水粒子大小;3)KDP,比差分相移,即单位距离上降水粒子导致的水平和垂直方向回波的相位差,主要反映了液态含水量。
2.Z-R关系:雷达反射率和降水之间的经验性关系,通常表述为,其中R为降水量,Z为雷达反射率,和为经验性参数,通常在不同地区及不同降水类型下有差异。
整体求解过程概述(摘要)
针对问题一,为了实现基于前一小时(10 帧)的实测雷达观测量(ZH、ZDR、KDP),对后续一小时(10 帧)的 ZH进行预报,本文首先建立了线性拟合与 RMSE 双驱动的局部离群因子检测(LOF)数据清洗模型,分别检测、去除了数据在时间和空间分布内的异常点,在数据清洗模型中,LOF 异常阈值设置为 LOF > 2。然后,本文系统性的分析了降雨各微物理参数对降雨气象模型、降雨散射参数、降雨雷达观测量产生影响的物理机理,说明了降雨气象雷达观测量随降雨量产生变化的本质是各尺寸雨滴粒子散射特性随不同雨滴谱相互耦合的结果。最后,基于降雨雷达观测量参数产生的物理机理,选用 NARX-RNN 模型构建神经网络,完成了基于前时间对后续时间雷达观测量的预估,结果表明相比较于测量真值,预估值相关性好(> 0.9),误差低(FRMSE < 1.7),说明该网络实现了对后续时间雷达观测量的预估,模型性能优异。
针对问题二,为了解决问题一模型中模型预报结果不清晰、不完备的模糊效应问题,本文首先对模糊效应产生机理进行了研究,分别对网络模型进行输入层信息挖掘和输出层信息细化,使网络输入层特征信息增加 336.67%,网络输出层信息增加 133.33%,实现了信息补足,缓解了模糊效应。在输入层上,本文引入雷达观测量三维梯度、降雨雷达强点、双偏振雷达类型参数作为模型补充输入,实现了输入层信息挖掘。在输出层上,本文分别从雷达回波细节和气象微物理参数两方面进行了扩充。在雷达回波细节扩充方面,本文基于双偏振雷达运作机理,形成 ZH,ZV,ZDR,KDP 双极化四输出完备雷达参数数据集,该数据集完整阐述了两个极化方向上的回波数据以及他们之间的内在关系。在气象微物理细节扩充方面,基于问题一降雨雷达观测量生成机理,推演计算生成包含气象类型、降水粒子大小、含水量信息的降雨气象微物理参数集,为含水量预估提供进一步支撑。最后,对扩充后的网络进行训练,完成了对上述输出层参数的预估,并与问题一进行了对比,结果表明网络预估相似度提升 5.37%,网络误差降低 30.54%,泛化性能提升 46.94%,网络性能显著提升。
针对问题三,为了实现基于雷达观测量(ZH、ZDR;不可使用 KDP),对降雨量实现定量估计,本文建立了测算融合的 BP 神经网络降水量反演模型。首先,基于问题一中讨论分析的降雨气象雷达观测量生成物理机理,利用自研软件,仿真生成不同降雨量下的气象散射特性,形成不同降雨量下的 ZH、ZDR仿真数据集,与测试真值对比结果表示仿真数据生成可靠(CC > 0.8; FRMSE < 2.0)。随后,将仿真与实测数据融合,进一步基于第二问,以雷达观测量扩充降雨气象微物理参数,形成测算融合数据集。最后,基于测算融合数据集,构建 BP 神经网络。结果表明,对降雨量的预估可靠性高(CC > 0.8; FRMSE < 2.0),网络性能优异。
针对问题四,为了对双偏振雷达在强对流降水中的贡献进行定量评估,并优化网络融合策略,实现模型针对降雨“突发性”、“局地性”的应对能力提升,本文分别从熵权-TOPSIS2评估模型构建和多模型网络融合两方面完成了问题要求。首先,本文基于第二问、第三问研究成果,形成 NARX-RNN-BP 多模型融合神经网络模型,作为后续评估基准网络模型。然后,调整基准模型输入层信息,分别对单偏振信息(仅 ZH)和不同等级的双偏振信息(仅ZDR、仅 KDP、ZH + ZDR和 ZH、ZDR、KDP 全参数)五种网络输入进行网络搭建,并对预估结果进行相关系数(CC)、均方根误差(RMSE)、相对偏置(RB)和分数均方差(FRMSE)四性能计算。随后,基于熵权-TOPSIS 模型构建评估模型,对五种网络输出进行评分,结果表明全参数相对评分 0.3520,排名第一,相比最低评分仅 ZH组提高 97.13%。接着,分别从局地性和突发性对网络进行优化。突发性方面,基于仿真数据构建了“是否有突发降雨”的随机森林分类模型,分类模型准确度达到 99.6%,完成了数据集的突发性标签制作。局地性方面,基于第二问降雨强点提取,对降雨强点构建 K-Means 分类模型,生成了局地雨区,完成了数据集的局地性标签制作。最后,对优化网络进行性能测试和评估,结果表明模型预估相似度极高(CC > 98%),误差极低(FRMSE < 1.2),相对模型分数提升 10.71% 到39.26% 不等,说明模型针对局地性与突发性泛化能力优异。
模型假设:
为了对进行准确、合理的描述,我们对关于降水临近预报模型做出以下假设:
(1)单位时间内可用于预测的数据量异常值数据不高于 10%,输入量在剔除异常值后仍满足网络最低需求;
(2)雷达系统在探测工作中的大气衰减造成的影响可以忽略;
(3)数据所使用的双偏振雷达的性能参数与某个“标准单偏振雷达”ZH值相同。
(4)雷达观测数据不因昼夜、气温因素造成偏差,并且探测范围内不存在非气象干扰(如鸟群、飞机等);
(5)雷达观测数据与降水实况之间存在稳定的统计关系;
(6)雷达系统的探测范围内仅存在水粒子气象;
(7)在不同的降雨中,相同尺寸的雨滴散射特性与雷达参量相同。
问题分析:
问题一需要根据附件一文件夹“NJU_CPOL_update2308”,下属包含“DBZ”,“KDP”,“ZDR”等次级文件,基于附件内数据建立可提取用于强对流临近预报双偏振雷达资料中微物理特征信息的数学模型。实现基于前面一小时(10 帧)的雷达观测量(ZH、ZDR、KDP),对后续一小时(10 帧)的 ZH 进行预报。观察题目所给附件分析,附件内包含了不同雷达作用距离下的雷达针对不同次降雨观测的随时间帧变化的网格化雷达实测数据。因为雷达实测数据受各种因素影响,在进行数学建模分析预测前首先要对数据进行清洗。随后,由于数据量巨大充足,可以根据“DBZ”,“KDP”,“ZDR”三数据随时间帧的变化规律基于神经网络实现未来雷达参数预测。
问题二指出,在常规模型构建中,往往会存在“倾向回归于平均”和“模糊效应”的问题,使得模型不清晰,输入输出不完备。例如在题目一中,仅将前十帧雷达观测量(ZH、ZDR、KDP)作为了网络输入,而网络输出仅有的后续一小时(10 帧)的 ZH值,这时的结果模糊、不完善。因此,本文理解问题二需要在问题一模型的基础上,基于附件一、二文件夹“NJU_CPOL_update2308”和“NJU_CPOL_kdpRain”,进一步挖掘数据,扩充模型输入,使网络模型输入更有条理,更有深度。同时,进一步挖掘双偏振雷达参数与参数之间,参数与降雨微物理特性之间的物理机理,使模型预报的结果在雷达回波纬度和气象微物理参数纬度两方面都更清晰,更细节。最终形成输入挖掘-输出细化的双驱动去模糊化临近预报模型。
问题三需要基于附件一、二文件夹“NJU_CPOL_update2308”中下属包含的“DBZ”“ZDR”两组次级文件和“NJU_CPOL_kdpRain”文件,建立可进行定量降水估计模型。实现基于的雷达观测量(ZH、ZDR;不可使用 KDP),对后降水量进行定量估计。在第二问中,我们已经探讨明白,只用 ZH、ZDR 进行降雨量预测会产生“模糊效应”,因此还需要对模型进行扩展。因此,本文理解首先基于物理机理研究,利用自研的复杂背景环境散射特性仿真软件,生成不同降雨量下的气象散射特性,并进一步形成不同降雨量下的 ZH、ZDR 仿真数据集。然后,将生成的仿真数据集与题目提供的测量数据集,在第二问和仿真软件的基础上,聚焦 ZH、ZDR 对气象微物理参数的影响,最终建立基于测算融合的 BP 神经网络降水量反演模型。仿真数据生成:基于降雨气象模型和降雨散射特性模型,利用自研的复杂背景环境散射特性仿真软件,仿真生成不同降雨量下的气象散射特性,形成不同降雨量下的 ZH、ZDR仿真数据集,为后续测算融合提供数据基础。测算融合 BP 神经网络:基于仿真获取的不同降雨量下的 ZH、ZDR 仿真数据集,融合题目附件提供的降雨实测数据,形成 ZH、ZDR测算融合数据集。然后在第二问和仿真软件对物理机理的探讨基础上,聚焦 ZH、ZDR 对气象微物理参数的影响,进一步构建 BP 神经网络,实现基于雷达观测量(ZH、ZDR)的降雨量分析模型。
问题四需要设计一种数学模型,使得可以对双偏振雷达参数在强对流降水预报中的贡献进行定量表述和比较。然后,基于这一种评价方式,优化模型数据的融合策略,实现模型对“突发性降雨类型”和“地域信息”降雨的泛化能力提升。因此,本文首先分别在考虑单偏振信息(仅 ZH)和考虑双偏振信息(ZH、ZDR 和 KDP)的条件下,融合 2、3 问的网络模型,分别对未来时间降雨量进行临近预报。然后,基于熵权-TOPSIS 方法,对相关系数(CC)、均方根误差(RMSE)、相对偏置(RB)和分数均方差(FRMSE)四种网络系数进行权重赋值,得到对预报结果的评估模型,评估双偏振雷达在强对流临近降雨预报中的贡献。最后,基于仿真数据和随机森林,建立突发性降雨分类模型,对附件数据库进行分类标签。然后,在神经网络中输入层加入考虑不同位置、降雨突发与否参数,实现对突发性和局地性强的强对流天气的预估。评估模型建立与双偏振雷达贡献度评估:分别在考虑单偏振信息(仅 ZH、仅 ZDR 和仅KDP)和考虑双偏振信息(ZH、ZDR 和 KDP)的条件下,融合 2、3 问的网络模型,实现多模融合,建立基准模型。然后,分别导入单偏振/双偏振雷达信息,对未来时间降雨量进行临近预报。利用熵权-TOPSIS 方法,对相关网络的相关系数(CC)、均方根误差(RMSE)、相对偏置(RB)和分数均方差(FRMSE)四种网络系数进行权重赋值,实现对模型预估性能的评估。然后通过前面建立的评估模型,评估双偏振雷达在强对流临近降雨预报中的贡献。优化数据融合的高泛化性网络模型:基于仿真数据库和随机森林分类模型,建立突发性降雨分类模型,对附件数据库进行分类标签。然后,在神经网络中输入层加入考虑不同位置、降雨突发与否,实现对突发性和局地性强的强对流天气的预估。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码:
部分程序代码如下:
1. void boxingtu::on_calculate_clicked(){
2. tmpvctor_up.clear();
3. tmpvctor_down.clear();
4. dnow.clear();
5. rcs_d.clear();
6. power_d.clear();
7. double G1 = Gmax; // 设置 G1 的值为 Gmax
8. phi0 = 0.;// 设置 phi0 的初始值为 0
9. int jihua = beam_jihua; // 设置 jihua 的值为 beam_jihua
10. ui->progressBar->setVisible(true); // 设置进度条可见
11. ui->label_4->setVisible(true); // 设置 label_4 可见
12. ui->progressBar->setRange(0,jmax1 * 1.2); // 设置进度条范围
13. ui->progressBar->setWindowModality(Qt::WindowModal); // 设置进度条模态窗口
14. ui->progressBar->setValue(jmax1 * 0.1); // 设置进度条的初始值
15. // 开始计算
16. if(ywy_class==1){
17. if(cloud_class==1||cloud_class==2){
18. iCloud_RCS(f0, pt, G1, beam, pw, 30.); // 根据条件调用 iCloud_RCS 函数
19. ui->progressBar->setValue(jmax1 * 0.2); // 更新进度条的值
20. }else{
21. wCloud_RCS(f0,pt,G1 ,beam,pw,30.); // 根据条件调用 wCloud_RCS 函数
22. ui->progressBar->setValue(jmax1 * 0.2); // 更新进度条的值
23. }
24. }else if(ywy_class==2){
25. fog_RCS(f0,pt, G1 ,beam,pw,30.); // 根据条件调用 fog_RCS 函数
26. ui->progressBar->setValue(jmax1 * 0.2); // 更新进度条的值
27. }else{
28. Rain_RCS(f0,pt,G1 ,beam,pw,jihua,30.); // 根据条件调用 Rain_RCS 函数
29. ui->progressBar->setValue(jmax1 * 0.2); // 更新进度条的值
30. }
31. dmax.clear(); // 清空 dmax
32. distancemax = 0; // 将 distancemax 的值设为 0
33. int jmax = jmax1; // 将 jmax 的值设为 jmax1
34. double j_delt = 1./ quick_f; // 将 j_delt 的值设为 1/quick_f
35. jmax1 << endl; // 输出 jmax1 的值并换行
36. missile m1 = tmpm; // 将 m1 的值设为 tmpm
37. double Tt = 0; // 将 Tt 的值设为 0
38. double dmax_ls = 0; // 将 dmax_ls 的值设为 0
39. ofstream rcsd; // 创建一个 ofstream 对象 rcsd
40. rcsd.open("rcsd.txt", ios::trunc); // 打开 rcsd.txt 文件,使用 trunc 模式
41. ofstream powerd; // 创建一个 ofstream 对象 powerd
42. powerd.open("powerd.txt", ios::trunc); // 打开 powerd.txt 文件,使用 trunc 模式
43. m1.cor = { m1.cor.x + Tt * m1.v.x, m1.cor.y + Tt * m1.v.y, m1.cor.z + Tt * m1.v.z };
// 更新 m1 的 cor 的值
44. if (ywy_class == 1){
45. cali = 1./3. * pi * pow((30. * sin(beam * pi/180.)),2) * 30. * cos(beam * pi/180.)
* pow(fsxs,3); // 根据条件计算 cali 的值
46. }
47. else{
48. cali = 1./3. * pi * pow((30. * sin(beam * pi/180.)),2) * 30. * cos(beam * pi/180.)
* pow(fsxs,3); // 根据条件计算 cali 的值
49. }
50. }
51.
52. for (int i = ceil(output_vector[rows - 1][0] / timestp); i < int(pulse_time / timestp);
i++) {
53. tmpvctor_down.push_back({Tt + j * j_delt + timestp * (i + 1),0.0}); // 将时间和
0.0 加入 tmpvctor_down 向量中
54. }
55. dmax.push_back(dmax_ls); // 将 dmax_ls 加入 dmax 向量中
56. srand(rand()); // 设置随机数种子
57. double u1 = (rand() % 100) * 0.01; // 生成 0 到 1 之间的随机数
58. double x_ave = 0.05 * u1 + 1 - 0.05; // 生成均匀分布的粒子
59. x_ave = x_ave * kkk; // 将粒子与 kkk 相乘
60. rcs_d.push_back(10 * log10(ans1.rcs * x_ave)); // 将 ans1.rcs 和 x_ave 相乘,然后取对
数并乘以 10,再加入 rcs_d 向量中
61. power_d.push_back(ans1.p * x_ave); // 将 ans1.p 和 x_ave 相乘,再加入 power_d 向量中
62. dnow.push_back(sqrt(move_x * move_x + move_y * move_y + move_z * move_z)); // 计算
move_x、move_y 和 move_z 的平方和的平方根,再加入 dnow 向量中
63. rcsd << sqrt(move_x * move_x + move_y * move_y + move_z * move_z) << " " << (10 *
log10(ans1.rcs * x_ave)) << endl; // 将 dnow 和 rcs_d 的最后一个元素输出到 rcsd 文件中
64. powerd << sqrt(move_x * move_x + move_y * move_y + move_z * move_z) << " " << ans1.p*
x_ave << endl; // 将 dnow 和 power_d 的最后一个元素输出到 powerd 文件中
65. if (distancemax < sqrt(move_x * move_x + move_y * move_y + move_z * move_z)){
66. distancemax = sqrt(move_x * move_x + move_y * move_y + move_z * move_z); // 如
果 distancemax 小于平方和的平方根,将其更新为平方和的平方根
67. }
68. double rcsnow = 10 * log10(ans1.rcs * x_ave), pnow = ans1.p * x_ave; // 计算 rcsnow
和 pnow
69. if (n == 0){
70. rcsmax = rcsnow; // 如果 n 等于 0,将 rcsmax 更新为 rcsnow
71. rcsmin = rcsnow; // 如果 n 等于 0,将 rcsmin 更新为 rcsnow
72. pmax = pnow; // 如果 n 等于 0,将 pmax 更新为 pnow
73. pmin = pnow; // 如果 n 等于 0,将 pmin 更新为 pnow
74. }
75. else{
76. if(rcsmax < rcsnow){
77. rcsmax = rcsnow; // 如果 rcsmax 小于 rcsnow,将其更新为 rcsnow
78. }
79. if(pmax < pnow){
80. pmax = pnow; // 如果 pmax 小于 pnow,将其更新为 pnow
81. }
82. if(rcsmin > rcsnow){
83. rcsmin = rcsnow; // 如果 rcsmin 大于 rcsnow,将其更新为 rcsnow
84. }
85. if(pmin > pnow){
86. pmin = pnow; // 如果 pmin 大于 pnow,将其更新为 pnow
87. }
1. void mkconebackground(vector<vector<double>>& yuandwusanweizuobiao) {
2. srand(zz); // 设置随机数种子
3. double totsize = maxcd; // 总大小
4. double length = pow(totsize, 3); // 总体积
5. rnum = RAND_MAX; // 随机数最大值
6. double aaa;
7. int h = 1;
8. double conemax = abs(tmpm.v.x);
9. if (abs(tmpm.v.y) > conemax) {
10. conemax = abs(tmpm.v.y);
11. }
12. if (abs(tmpm.v.z) > conemax) {
13. conemax = abs(tmpm.v.z);
14. }
15. double vcone[3] = { tmpm.v.x / conemax, tmpm.v.y / conemax, tmpm.v.z /
conemax }; // 轴线向量
16. double theta = beam * pi / 180.;
17. double pc0[3] = { vcone[0] + pos0_x,vcone[1] + pos0_y ,vcone[2] +
pos0_z };//
18. ofstream coneshow;
19. coneshow.open("coneshow.txt", ios::trunc);
20. double kbei;
21. if (ywy_class == 3) {
22. kbei = 5000. / (10 * log10(rain1.R * 100 + 1));
23. }
24. else if (ywy_class == 2) {
25. kbei = (100. * log10(fog1.V * 1000 + 1));
26. }
27. else {
28. kbei = 10.;
29. }
30. for (int i = 0; i < ceil((length * 0.2e6 / kbei) / (maxcd * maxcd * maxcd));
i++) { // 循环空间粒子总数次
31. aaa = rand();
32. double x = (aaa / rnum) * totsize;
33. aaa = rand();
34. double y = (aaa / rnum) * totsize;
35. aaa = rand();
36. double z = (aaa / rnum) * totsize;
37. double pc1[3] = { x,y,z };
38. coneshow << x << " " << y << " " << z << endl;
39. yuandwusanweizuobiao.push_back({ x,y,z });
40. if (ywy_class == 3) {
41. for (int iii = 1; iii < 50; iii++) {
42. yuandwusanweizuobiao.push_back({ x,y + 0.05 * iii,z });
43. coneshow << x << " " << y + 0.05 * iii << " " << z << endl;
44. }
45. }
46. }
47. cali = h * kbei;
48. coneshow.close();
49. cout << "moxjianli _ ed" << endl;
50. rangexmax = totsize; rangeymax = totsize; rangezmax = totsize;
51. rangexmin = 0; rangeymin = 0; rangezmin = 0;
52. }