一种基于物理损失引导神经网络的湖冰厚度预测方法

本发明涉及湖泊冰情遥感与预测，具体为一种基于物理损失引导神经网络的湖冰厚度预测方法。

背景技术：

1、浅水湖泊冰厚是冬季湖泊热量交换、冰情评估和冰季物候分析的重要状态变量。受观测成本、冰面安全和站点分布限制，连续湖冰厚度观测通常较为稀疏。传统过程模型能够表达能量收支和冰雪热力过程，但需要较多参数及站点化配置；纯数据驱动模型能够利用气象变量进行预测，却容易在小样本和季节滚动预测中产生不符合物理规律的波动。因此，需要一种能够处理不等间隔时间序列、利用多湖泊样本并在训练过程中引入多时间尺度物理约束的湖冰厚度预测方法。

技术实现思路

1、本发明的目的在于提供一种基于物理损失引导神经网络的湖冰厚度预测方法，以解决上述背景技术中提出的问题。

2、为实现上述目的，本发明提供如下技术方案：

3、一种基于物理损失引导神经网络的湖冰厚度预测方法，包括以下步骤：

4、s1、数据底座构建：获取目标湖泊及辅助训练湖泊的数据；对所述数据进行时间戳标准化、时间一致性处理、变量单位统一和缺失值处理，并生成前一时刻冰厚状态变量、度日变量、湖泊深度变量和雪深代理变量；

5、s2、训练样本与预测窗口构建：形成带有相对时间通道的多变量时间序列窗口；

6、s3、连续时间控制路径生成：将所述多变量时间序列窗口转换为连续时间控制路径的插值系数；所述连续时间控制路径采用hermite三次插值表示；

7、s4、神经控制微分方程编码与非负冰厚输出：将所述连续时间控制路径输入神经控制微分方程模型，并将最终潜在状态映射为非负湖冰总厚度预测值；

8、s5、季节尺度度日物理损失构建：构建季节尺度参考冰厚曲线，用于约束湖冰形成、增长和消融的季节演化包络；

9、s6、日尺度能量-冰厚物理损失构建：计算表面能量收支、水-冰界面热通量、能量限制增长增量和能量限制消融增量，用于约束每日冰厚变化的方向和幅度；

10、s7、复合物理损失函数训练：将实测冰厚拟合损失、季节尺度度日物理损失和日尺度能量-冰厚物理损失组合为复合训练目标，对所述神经控制微分方程模型进行物理损失引导训练；

11、s8、自回归季节滚动预测：在完成训练后，使用目标湖泊逐日气象强迫数据进行季节滚动预测，将上一日预测得到的湖冰厚度反馈为下一日预测的前一时刻冰厚状态变量，并结合开放水面判定、能量预算一致性判定和日增长合理性判定，生成连续日尺度湖冰厚度轨迹；

12、s9、湖冰物候提取与模型评价：计算冰厚预测精度指标、季节轨迹一致性指标和物候误差指标；

13、s10、结果输出。

14、进一步的，该湖冰厚度预测方法，包括以下步骤：

15、s1、数据底座构建：以目标湖泊为预测对象，获取目标湖泊及辅助训练湖泊的实测湖冰厚度数据、湖泊深度先验数据和气象强迫数据；对所述数据进行时间戳标准化、时间一致性处理、变量单位统一和缺失值处理，并生成前一时刻冰厚状态变量、度日变量、湖泊深度变量和雪深代理变量；所述雪深代理变量按式（1）生成：

16、（1）

17、其中，为第日雪深代理变量，为雪深代理变量上限，为水密度，为雪密度，为滚动累计窗口天数，为第日降雪水当量；当求和窗口中的降雪量为负值时取零；

18、s2、训练样本与预测窗口构建：针对每一个训练样本日期或待预测日期，按照预设回看窗口同步提取气象强迫变量和前一时刻冰厚状态变量，形成带有相对时间通道的多变量时间序列窗口；

19、s3、连续时间控制路径生成：将所述多变量时间序列窗口转换为连续时间控制路径的插值系数，使模型能够沿连续时间路径读取不等间隔或稀疏观测条件下的气象-冰厚状态历史；所述连续时间控制路径采用hermite三次插值表示；

20、s4、神经控制微分方程编码与非负冰厚输出：将所述连续时间控制路径输入神经控制微分方程模型，由所述神经控制微分方程模型沿连续时间控制路径更新潜在状态，并将最终潜在状态映射为非负湖冰总厚度预测值；当表示时间区间]上的连续时间控制路径时，所述神经控制微分方程模型按式（2）至式（4）进行编码和输出：

21、（2）

22、（3）

23、（4）

24、其中，为潜在状态，为初始状态映射网络，为向量场网络，和为输出读出层参数，为第个样本的非负湖冰总厚度预测值；

25、所述向量场网络由多层感知机参数化，所述多层感知机包括线性层、relu激活层和dropout层；所述神经控制微分方程采用rk4求解器积分，目标变量为以米为单位的湖冰总厚度且不进行对数变换；

26、s5、季节尺度度日物理损失构建：基于近地面空气温度、冰季起始月份和冻结温度阈值，计算累积冻结度日afdd和累积融化度日atdd，并构建季节尺度参考冰厚曲线，用于约束湖冰形成、增长和消融的季节演化包络；

27、s6、日尺度能量-冰厚物理损失构建：基于前一时刻冰厚、相邻观测时间间隔、气象强迫数据、湖泊深度先验数据和雪深代理变量，计算表面能量收支、水-冰界面热通量、能量限制增长增量和能量限制消融增量，用于约束每日冰厚变化的方向和幅度；

28、s7、复合物理损失函数训练：将实测冰厚拟合损失、季节尺度度日物理损失和日尺度能量-冰厚物理损失组合为复合训练目标，对所述神经控制微分方程模型进行物理损失引导训练，使模型同时满足观测拟合、季节演化和日尺度能量一致性；

29、s8、自回归季节滚动预测：在完成训练后，使用目标湖泊逐日气象强迫数据进行季节滚动预测，将上一日预测得到的湖冰厚度反馈为下一日预测的前一时刻冰厚状态变量，并结合开放水面判定、能量预算一致性判定和日增长合理性判定，生成连续日尺度湖冰厚度轨迹；对于季节首个预测日，在无可用历史冰厚观测时，将前一时刻冰厚初始化为0，将前一时刻冰厚状态可用性标识设为无效，并将首日相邻时间间隔设为0；按照预设稀疏窗口回退规则构造首日输入窗口；当存在指定历史观测时，优先采用该历史观测作为前一时刻冰厚状态。

30、s9、湖冰物候提取与模型评价：基于所述连续日尺度湖冰厚度轨迹和预设冰存在阈值提取冻结开始、完全消融和冰盖持续时间中的一种或多种湖冰物候指标，并计算冰厚预测精度指标、季节轨迹一致性指标和物候误差指标；

31、s10、结果输出：输出训练后的神经控制微分方程模型、日尺度湖冰厚度预测序列、湖冰物候指标、模型评价指标、训练配置、物理损失诊断信息和预测日志。

32、进一步的，所述气象强迫数据包括10m高度风速、2m空气温度、相对湿度、下行短波辐射、下行长波辐射、海平面气压、地表气压、降水量和降雪量中的一种或多种；所述变量单位统一包括将辐射累积量换算为瓦特每平方米，将降水量和降雪量换算为日尺度水当量，并将负的降水量或降雪量裁剪为零。

33、进一步的，所述前一时刻冰厚状态变量包括前一时刻湖冰厚度、距前一时刻冰厚状态的时间间隔和前一时刻冰厚状态可用性标识；所述预设回看窗口为7日回看窗口，窗口至少包含两个时间点，并在训练阶段锚定于实测湖冰厚度日期，在季节滚动预测阶段锚定于逐日预测日期。

34、进一步的，所述实测冰厚拟合损失采用式（5）至式（7）的huber损失：

35、（5）

36、（6）

37、（7）

38、其中，为第个样本的实测湖冰厚度，为huber损失阈值，为样本数。预测湖冰厚度以表示，实测湖冰厚度以表示。

39、进一步的，所述累积冻结度日afdd和累积融化度日atdd按湖泊和冰季分别累计，并按式（8）至式（10）计算：

40、（8）

41、（9）

42、（10）

43、其中，为冻结温度阈值，为第日近地面空气温度，为第日所属冰季内截至第日的样本集合；当相邻记录跨越不同冰季时，对应时间间隔不参与上一冰季累计；其中，为以天计的相邻记录时间间隔，为以秒计的相邻记录时间间隔，且=86400；度日累计使用，能量增量式使用。

44、所述季节尺度参考冰厚曲线按式（11）计算：

45、（11）

46、其中，为第日日度参考冰厚，为冻结增长系数，为融化衰减系数；和为正值参数，并通过softplus变换或其他正值约束方式保持非负；

47、所述季节尺度度日物理损失按式（12）计算：

48、（12）

49、其中，为增长阶段且具有稳定冰层和可用度日驱动的样本集合，为消融阶段且具有稳定冰层和可用度日驱动的样本集合，和分别为增长阶段和消融阶段度日损失权重。

50、进一步的，所述日尺度能量-冰厚物理损失中，表面反照率、饱和水汽压、比湿和表面能量项按式（13）至式（21）计算：

51、（13）

52、（14）

53、（15）

54、（16）

55、（17）

56、（18）

57、（19）

58、（20）

59、（21）

60、其中，为表面反照率，为雪面反照率，为冰面反照率，为参考雪深，为下行长波辐射，为有效短波辐射，为冰面发射率，为stefan-boltzmann常数，为空气密度，为空气定压比热，和分别为感热和潜热交换系数，为风速，为升华潜热，为降水热项；其中，空气水汽压=()，表面水汽压=()；压力以pa为单位，并满足−0.378>0；参考雪深为正值，或在计算时采用max(，ε)，其中ε为预设正数。

61、所述有效短波辐射通过本地正午短波辐射与日平均修正因子相乘得到，并按式（22）至式（25）计算：

62、（22）

63、（23）

64、（24）

65、（25）

66、其中，为日序号，为湖泊纬度对应的弧度值，为太阳赤纬，为日落时角，和分别为日平均修正因子的下限和上限；当预设参数直接给定时，采用所述预设参数作为日平均修正因子。

67、进一步的，所述水-冰界面热通量按式（26）至式（30）计算：

68、（26）

69、（27）

70、（28）

71、（29）

72、（30）

73、其中，为冰下有效水深，为湖泊深度先验，为前一时刻湖冰厚度，和分别为冰下有效水深下限和上限，、和分别为雪、冰和水体的短波衰减参数，为进入冰雪层并向水体透射的短波比例，和为水-冰界面热通量参数，为向上水-冰界面热通量的非负部分；

74、所述能量限制增长增量和能量限制消融增量按式（31）至式（39）计算：

75、（31）

76、（32）

77、（33）

78、（34）

79、（35）

80、（36）

81、（37）

82、（38）

83、（39）

84、其中，为相邻冰厚状态的时间间隔，为冰密度，为冰融化潜热，、、和分别为增长、表面消融、底部消融和升华相关质量损失效率系数，为雪保温系；

85、所述日尺度能量-冰厚物理损失按式（40）至式（42）计算：

86、（40）

87、（41）

88、（42）

89、其中，为满足局部能量约束启用条件且大于零的样本集合，为满足局部能量约束启用条件且不大于零的样本集合，和分别为增长和消融松弛阈值，和分别为方向性惩罚权重，和分别为增长约束和消融约束权重；

90、所述局部能量约束启用条件包括：前一时刻冰厚状态可用，湖泊深度先验为有限值，所述相邻冰厚状态的时间间隔大于零且不超过预设最大间隔，度日驱动数据可用，并且度日参考冰厚大于零；

91、所述复合训练目标按式（43）计算：

92、（43）

93、其中，为实测冰厚拟合损失权重，为实测冰厚拟合损失，为季节尺度度日物理损失，为日尺度能量-冰厚物理损失；

94、所述自回归季节滚动预测中，对于第个预测日，将第个预测日的湖冰厚度预测值作为前一时刻湖冰厚度，将相邻预测日间隔作为前一时刻冰厚状态的时间间隔，将前一时刻冰厚状态可用性标识设为有效；当开放水面判定满足气温高于预设阈值且前一日预测冰厚低于预设阈值时，将预测冰厚置于开放水面状态或进行非负裁剪；当能量预算一致性判定显示预测日增长或消融量超过预设能量限制容差时，对预测冰厚变化幅度进行限制；对于季节首个预测日，在无可用历史冰厚观测时，将前一时刻冰厚初始化为0，将前一时刻冰厚状态可用性标识设为无效，并将首日相邻时间间隔设为0；按照预设稀疏窗口回退规则构造首日输入窗口；当存在指定历史观测时，优先采用该历史观测作为前一时刻冰厚状态。

95、所述湖冰物候指标按式（44）和式（45）提取：

96、（44）

97、（45）

98、其中，为预设冰存在阈值，为有冰日标识；冻结开始为冰季内首次满足持续有冰条件的日期；完全消融为冰季末最后一个有冰日之后的首个无冰日期；冰盖持续时间为与之间的日数差。

99、进一步的，还包括消融实验评价步骤：在保持输入变量、回看窗口、插值方式、神经控制微分方程结构和训练协议一致的条件下，分别训练完整物理损失模型、仅实测冰厚拟合模型、保留季节尺度度日物理损失且移除日尺度能量-冰厚物理损失的模型、以及保留日尺度能量-冰厚物理损失且移除季节尺度度日物理损失的模型，并比较不同模型在冰厚预测精度、季节滚动轨迹和湖冰物候指标上的差异。

100、一种基于物理损失引导神经网络的湖冰厚度预测系统，包括处理器、存储器以及存储在所述存储器上并能够由所述处理器执行的程序模块，所述程序模块包括数据预处理模块、训练样本与预测窗口构建模块、连续时间控制路径生成模块、神经控制微分方程预测模块、季节尺度度日物理损失模块、日尺度能量-冰厚物理损失模块、复合损失训练模块、自回归季节滚动预测模块、湖冰物候提取与模型评价模块以及结果输出模块；所述处理器执行所述程序模块时，实现所述的方法。

101、一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，实现所述的方法。

102、与现有技术相比，本发明的有益效果是：

103、本发明通过连续时间控制路径对不等间隔的气象—冰厚历史进行编码，能够利用稀疏观测样本构建时序预测模型；通过非负读出层约束预测冰厚的取值范围；通过季节尺度度日物理损失约束冰季演化相位，并通过日尺度能量—冰厚物理损失约束相邻时刻冰厚变化的方向和幅度，以降低季节滚动预测中不符合物理过程的波动；能够输出连续冰季冰厚轨迹以及冻结开始、完全消融和冰盖持续时间等物候指标，并记录训练配置和物理损失诊断信息。