文章目录
- 介绍
- 代码
- 参考
介绍
冰川湖溃决洪水（GLOF）是山区面临的一大灾害，但关于其频率在近几十年内是否有所增加以及这种趋势与气候变化之间的关联程度，仍存在相当大的不确定性。在此，我们基于冰碛湖溃决洪水这一新视角，对全球 1900 年至 2020 年期间的 609 起事件进行了分析。历史报告和地貌学证据提供的信息表明，在 20 世纪 70 年代之前，全球报告的 GLOF 频率呈现低但波动上升的趋势。然而，自 20 世纪 80 年代以来，这一趋势显著加速，1981 年至 1990 年期间的年均 GLOF 数量为 5.2 起，而 2011 年至 2020 年期间则上升至 15.2 起。总体而言，报告的 GLOF 频率的长期趋势与全球气温的变化密切相关，呈现出大约 20 年的滞后相关模式。我们引入了 GLOF 响应时间的概念来解释这种延迟反应，这种反应归因于气候变暖导致的冰川退缩、冰川湖扩张以及此类湖泊周围的斜坡失稳，最终引发了 GLOF。
由冰碛湖堰塞形成的冰川湖溃决洪水（GLOF）是山区地区的一大安全隐患1。与之相关的连锁反应往往跨越数十至数百公里，并导致深远的社会和地貌后果2,3。例如，在 2023 年的锡金喜马拉雅地区，南隆纳克 GLOF 造成 55 人死亡，另有 74 人失踪，并在下游超过 160 公里的范围内造成了严重破坏，摧毁了约 30 座桥梁、2000 座建筑物和三个水电站4，5。全球范围内还有其他严重的 GLOF 现象已被详细记录，包括在加拿大不列颠哥伦比亚省、秘鲁科尔迪拉布拉纳山脉以及巴塔哥尼亚安第斯山脉等地，并且正受到越来越多的科学和政策关注10。GLOF 的发生通常与冰川退缩（即冰川萎缩导致山谷坡面崩塌以及不稳定冰前沿的冰块滑入湖泊）、永久冻土退化（岩崩、滑坡和冰碛坝中死亡冰的融化）以及极端天气事件（极端降雨和融水事件）有关，而所有这些都因全球变暖而加剧11，12，13。然而，全球性冰川湖溃决洪水与气候变化之间的关联仍存在争议，并且常常受到潜在遗漏和记录错误的影响14。对现有全球性冰川湖溃决洪水清单的分析迄今表明，自 20 世纪 80 年代以来，其发生频率在区域和全球范围内均出现了下降或停滞现象15、16、17，这与人们直观的预期相矛盾，即鉴于近期的强烈变暖18 以及冰川湖数量的增加19、20，全球性冰川湖溃决洪水的发生本应有所增加。
不断更新和修订的全球冰川湖溃决洪水（GLOF）清单对于解决这些差异至关重要。它们能够对 GLOF 的发生频率和与气候的关联进行分析，并为更广泛的应用提供信息，例如验证 GLOF 易发性评估的可靠性21、确定用于 GLOF 情景建模的湖泊22、23、24，以及为 GLOF 风险缓解和政策制定提供信息25、26。可以使用两种方法构建精确的 GLOF 时间序列。第一种方法是汇集文献来源中报告的事件，包括期刊论文、新闻报道和当地行政记录27。第二种方法结合了灾后地貌分析，利用 GLOF 的痕迹，如湖泊水位下降、堤坝溃决、冲积扇和下游破坏等，来识别以前未报告的事件28、29。对来自不同来源的图像进行手动或计算机自动分析可以极大地丰富清单。喜马拉雅地区17、热带安第斯山脉7 和南安第斯山脉30 的研究表明，仅依靠被动报告的信息会低估冰川湖溃决洪水的量 0.5 至 2 倍。因此，构建全球冰川湖溃决洪水（GLOF）数据库需要持续更新已报告的事件，并通过地貌评估主动识别未报告的历史性 GLOF。此外，对于模糊的情况，GLOF 的痕迹可以确认其可靠性31，而像 Landsat 和 Sentinel 这样的广泛卫星图像档案则可用于进一步限定那些爆发时间窗口较大的事件的 GLOF 时间。最新的全球数据库（版本 4.1）包含了 1900 年至 2020 年间发生的 463 个冰碛湖溃决洪水事件32，但仍需要在完整性、连续性和准确性方面进行改进。
在这项研究中，我们更新了关于冰碛坝型泥石流的全球统计信息，从而能够对它们的分布、发生频率以及与气候的关系进行系统分析。我们采用了由数据集整合、爆发时间校准以及未报告泥石流的识别组成的三步框架（见方法部分）来提高数据集的质量，并支持在不同空间和时间尺度上对泥石流特征进行可靠的评估。主要发现表明，自 20 世纪以来，泥石流的频率呈阶段性上升趋势，这主要由全球变暖加剧所驱动。
代码
https://github.com/Taigang1/GLOF/blob/main/GAM.py
importpandasaspdimportnumpyasnpimportmatplotlib.pyplotaspltfrompygamimportPoissonGAM,simportstatsmodels.apiassmimportstatsmodels.formula.apiassmffromscipyimportstats# 读取Excel文件file_path='H:\P20230801_Moraine-dammed lake outburst activities\GAM model\Trend.xlsx'# 替换为你的文件路径data=pd.read_excel(file_path)# 提取自变量和因变量year=data['Year'].values frequency=data['Frequency'].values# 定义候选的自由度范围df_candidates=range(3,15)# 可调控，很重要，确保 n_splines > spline_order# 存储每个自由度对应的AIC值和交叉验证误差aic_values=[]cv_errors=[]smooth_term_p_values=[]estimated_slopes=[]fordfindf_candidates:print(df)try:# 拟合GAM模型，spline_order是控制平滑项的多项式阶数的参数，lam控制平滑度的参数gam=PoissonGAM(s(0,n_splines=df,spline_order=2),lam=0.1)gam.gridsearch(year[:,None],frequency)# 计算AICaic_values.append(gam.statistics_['AIC'])# 交叉验证cv_error=gam.score(year[:,None],frequency)cv_errors.append(cv_error)# 获取光滑项的p值和斜率p_values=gam.statistics_['p_values']slopes=gam.coef_ smooth_term_p_values.append(p_values)estimated_slopes.append(slopes)exceptExceptionase:print(f"Error with df={df}:{e}")# 选择最优的自由度best_df_aic=df_candidates[np.argmin(aic_values)]best_df_cv=df_candidates[np.argmax(cv_errors)]print(f"Best df by AIC:{best_df_aic}")print(f"Best df by CV:{best_df_cv}")# 选择最优自由度进行最终模型拟合（例如，选择交叉验证结果）best_df=best_df_cv# 拟合最终GAM模型gam=PoissonGAM(s(0,n_splines=best_df,spline_order=2),lam=0.1)gam.gridsearch(year[:,None],frequency)# 预测值和置信区间year_pred=np.linspace(year.min(),year.max(),1000)# 增加样本点数拟合置信区间pred_mean=gam.predict(year_pred[:,None])# 计算置信区间confidence_intervals=gam.confidence_intervals(year_pred[:,None],width=0.95)# 提取下限和上限lower_bound=confidence_intervals[:,0]upper_bound=confidence_intervals[:,1]# 获取统计指标（估计自由度反映模型的复杂度和光滑程度。较高的edf表示模型更复杂）p_values=gam.statistics_['p_values']edf=gam.statistics_['edof']aic=gam.statistics_['AIC']deviance=gam.statistics_['deviance']#较低的AIC值表示模型在平衡拟合优度和复杂度方面表现较好。# 计算χ²统计量（较低的统计量表示拟合较好）observed=frequency predicted=gam.predict(year[:,None])chi2=np.sum((observed-predicted)**2/predicted)# 计算平均斜率（反映曲线整体趋势）average_slope=np.mean(gam.coef_)# 使用 statsmodels 进行线性回归来计算斜率的显著性# 将数据准备为适合 statsmodels 的格式df=pd.DataFrame({'year':year,'frequency':frequency})model=smf.poisson('frequency ~ year',data=df).fit()# 获取斜率的标准误差、t值和p值（判断斜率是否显著）avg_slope_std_error=model.bse['year']avg_slope_t_value=model.tvalues['year']avg_slope_p_value=model.pvalues['year']# 输出评估指标print(f"p-values:{p_values}")print(f"Estimated Degrees of Freedom (edf):{edf}")print(f"Chi-squared statistic (χ²):{chi2}")print(f"AIC:{aic}")print(f"Average Slope of Smooth Terms:{average_slope}")print(f"Average Slope Standard Error:{avg_slope_std_error}")print(f"Average Slope t-value:{avg_slope_t_value}")print(f"Average Slope p-value:{avg_slope_p_value}")# 输出光滑函数的p值和斜率print(f"Smooth term p-values:{p_values}")print(f"Estimated slopes of smooth terms:{gam.coef_}")# 将拟合曲线和不确定区间上下线的x和y坐标信息导出到Excel文件output_data=pd.DataFrame({'Year':year_pred,'Fitted':pred_mean,'Lower Bound (95%)':lower_bound,'Upper Bound (95%)':upper_bound})output_file_path='H:\P20230801_Moraine-dammed lake outburst activities\GAM model\GAM_Fit_Results.xlsx'# 替换为你想要的输出文件路径output_data.to_excel(output_file_path,index=False)# 将评估指标添加到Excel文件中withpd.ExcelWriter(output_file_path,mode='a',engine='openpyxl')aswriter:metrics_data=pd.DataFrame({'Metric':['p-values','Estimated Degrees of Freedom (edf)','Chi-squared statistic','AIC','Average Slope','Average Slope Standard Error','Average Slope t-value','Average Slope p-value'],'Value':[p_values,edf,chi2,aic,average_slope,avg_slope_std_error,avg_slope_t_value,avg_slope_p_value]})metrics_data.to_excel(writer,sheet_name='Metrics',index=False)print(f'Results saved to{output_file_path}')# 绘图plt.figure(figsize=(12,6))plt.plot(year,frequency,'o',label='Observed')plt.plot(year_pred,pred_mean,'r-',label='Fitted')plt.fill_between(year_pred,lower_bound,upper_bound,color='r',alpha=0.2,label='95% CI')plt.xlabel('Year')plt.ylabel('Frequency of Glacier Lake Outburst Floods')plt.title('GAM Fit with Optimal df for Frequency of Glacier Lake Outburst Floods')plt.legend()plt.show()参考
- High frequency of moraine-dammed lake outburst floods driven by global warming
- https://github.com/Taigang1/GLOF/blob/main/
