量化气候变化影响的不确定性：卡加延河流域的贝叶斯方法

《International Journal of Disaster Risk Reduction》：Quantifying uncertainty in climate change impacts: A Bayesian approach for the Cagayan River Basin

　　栗原勇太｜宫本守
日本东京国立政策研究研究生院（GRIPS），106-0032

**摘要**
气候变化会加剧连续多日的极端降雨，并增加破坏性的河流洪水。然而，罕见事件的降水具有高度不确定性，并以非线性的方式导致洪水和损失。本研究通过将贝叶斯极值推断与降雨-洪水-损失链条中的不确定性传播相结合，为菲律宾吕宋岛北部的卡加延河流域开发了一种概率性洪水影响评估方法。该研究使用贝叶斯广义极值推断对流域平均1天和7天的年最大降水量进行了建模。通过结合三个CMIP6全球气候模型（GCM），使用基于连续排序概率分数的贝叶斯模型平均法估计了SSP1–2.6和SSP3–7.0情景下的100年一遇降雨量（RL100）。洪水响应通过由固定风暴模板（Ondoy–Pepeng）驱动的降雨-径流-淹没（RRI）模型模拟来表示，并根据降雨总量进行缩放。从这些模拟中，得出了建筑物、水稻和玉米的降雨-峰值阶段关系以及降雨-损失函数。通过对RL100降雨量进行蒙特卡洛抽样，并通过这些关系传播，得到了峰值阶段和事件损失的后验分布，以中心（50%）和包含尾部（95%）区间的可信范围宽度进行了总结。SSP1–2.6情景下的预测结果与当前气候相当，不确定性相对较小，而SSP3–7.0情景下的预测结果将RL100向上移动并扩大了分布范围。当不确定性进一步转化为峰值阶段和损失时，建筑物的损失区间宽度几乎翻了一番（从19.44亿菲律宾比索增加到38.17亿菲律宾比索），水稻的损失区间宽度也有所增加。这些结果表明，基于确定性RL100的估计在更强的气候强迫下可能低估了尾部风险，因此概率性报告对于制定稳健的适应计划至关重要。

**1. 引言**
预计气候变化将加剧水文气象极端事件，增加高影响洪水的可能性，对基础设施、农业和生计产生重大影响[1]，[2]。在易发生洪水的地区，这为预先投资适应措施和风险降低提供了有力依据。然而，气候预测的不确定性很大，特别是对于导致严重河流洪水的罕见连续多日极端降雨，这使得决策变得复杂。在极端事件方面，广义极值（GEV）模型是年最大降水量的标准工具，而贝叶斯推断越来越多地被用来通过后验可信区间表示参数不确定性和尾部行为[3]，[4]。在气候模型层面，多模型之间的差异和模型间的独立性有限，这促使人们采用显式的基于技能的权重方法而不是简单的平均法[5]，[6]。使用诸如连续排序概率分数（CRPS）之类的适当评分规则进行的贝叶斯模型平均（BMA）为将预测结果合并为校准后的不确定性估计提供了一个连贯的框架[7]，[8]，[9]。将这种降雨不确定性转化为洪水损失需要通过基于物理的洪水模拟和损失函数进行传播。这些过程中的非线性效应可能会将尾部降雨量的微小差异放大为经济影响的巨大差异，使得基于确定性的设计值容易低估极端结果。

尽管有越来越多的研究，但从气候模型预测到极端降雨推断、基于物理的洪水模拟以及直接经济损失估计的端到端概率框架仍然较为罕见。许多研究仅涉及链条的一部分——降雨或洪水阶段——而没有扩展到经济影响，或者采用压制不确定性的确定性设计值。本研究针对这三个具体空白进行了探索：（i）气候预测中的多模型不确定性很少与基于GEV的参数不确定性结合成一个单一的后验分布；（ii）降雨在损失阶段的非线性放大效应（即尾部风险的增加）很少被明确量化；（iii）对于菲律宾最大且最易发生洪水的流域之一——卡加延河流域（CRB），在对比排放情景下的概率性洪水影响评估几乎不存在。本研究通过关注100年一遇多日极端降雨事件的不确定性及其对洪水阶段和直接损失的影响，填补了这些空白。使用贝叶斯广义极值推断对极端降雨进行了建模，结合了基于历史表现的CRPS衍生出的模型权重，并使用RRI模型和针对建筑物及主要作物的直接损失框架将不确定性传播到洪水响应[10]。[11]。该分析在链条的每个阶段都保留了概率分布，而不是简化为点估计，并报告了在不同气候情景下的降雨、洪水阶段和损失的可信范围。

**2. 研究区域和数据**
**2.1 研究区域**
研究区域是菲律宾吕宋岛北部的卡加延河谷地区（Region II）。卡加延河是菲律宾最长的河流，流域面积全国最大。年降雨量差异显著，北部低地的降雨量约为1000毫米，而南部山区则达到约3000毫米。

**2.2 水文气象观测和卫星降水数据**
降雨观测数据来自PAGASA站点，包括由河流洪水预报和预警中心（RFFWC）运营的长期日降水计和小时降水计；这些小时降水计中含有大量缺失值[12]。为了弥补覆盖范围的不足，使用了GSMaP v8 MVK卫星产品，并根据[12]进行了偏差校正。观测数据覆盖时间为1998–2024年（共27个年最大降水量）。

**2.3 气候模型数据（历史和情景）**
未来气候预测数据来自NEX-GDDP-CMIP6数据集中的CMIP6全球气候模型（分辨率为0.25°/天）。根据ERA5对上述降水及关键气象变量的评估（相关性和均方根误差[12]），选择了三个模型（ACCESS-CM2、CanESM5和EC-Earth3-Veg-LR）。在两个对比情景下进行了预测，以涵盖不确定性范围。历史参考期为1995–2014年，未来预测期为2040–2099年。为了校正日模拟降雨量的偏差，应用了基于分位数的三步统计降尺度方法，以偏差校正后的GSMaP降水数据作为参考（附录）。校正过程是在将GSMaP数据重新投影到模型网格后，并使用1998–2017年的数据进行。卫星记录始于1998年，而模型历史数据始于1995年；这种时间不匹配是由于数据可用性的限制，但假设重叠期足以表征当前气候的偏差。

**2.4 暴露和脆弱性数据**
暴露和脆弱性数据遵循为台风Ulysses开发的直接损失框架[11]。水稻和玉米的分布范围基于随机森林土地利用和土地覆盖分类[11]，建筑物的暴露范围则使用Open Buildings V3多边形表示[13]。住宅建筑[14]、[15]以及水稻和玉米[16]、[17]的深度-损失关系已确定。资产价值数据来自官方统计[18]，[19]。

**3. 方法**
**3.1 概述：概率性洪水风险和损失框架**
该框架通过将设计降雨量的不确定性通过洪水模拟和直接损失估计进行传播，来量化未来气候情景下100年一遇洪水事件的不确定性。首先，使用贝叶斯GEV框架对流域平均1天和7天的年最大降水量进行了建模。选择7天累积窗口作为主要风险指标，因为CRB地区最具破坏性的历史事件——2009年9月至10月的台风Ondoy和Pepeng——是由大约440毫米的7天累积降水量驱动的，远超过1天的最大降水量约140毫米。这证实了在该流域，持续累积而非短暂强度是极端河流洪水的主要原因。保留1天序列用于比较和描述性参考。其次，使用贝叶斯模型平均（BMA）表示气候模型的不确定性，模型权重基于历史CRPS表现得出。然后，使用固定观测风暴模板（Ondoy–Pepeng）对降雨-径流-淹没（RRI）模型进行驱动，并改变流域平均7天降水量。这些模拟用于推导降雨-峰值阶段和降雨-损失关系。最后，通过对RL100降雨量的蒙特卡洛抽样，并通过这些关系传播，得到峰值阶段和事件损失的后验分布。

**3.2 观测到的年最大降水量和贝叶斯设计降雨量**
使用移动窗口提取法为流域平均1天和7天累积降水量构建了年最大降水量序列。对于观测数据，序列覆盖时间为1998–2024年（共27个年最大降水量）；对于模型数据，分别为历史和未来时期提取了独立的序列。年最大降水量使用广义极值分布进行建模，即（1）F(x;μ,σ,ξ)=exp^(-1+ξ)(x-μ)/σ^(-1/ξ)，其中1+ξ((x-μ)/σ)>0，μ是位置参数，σ>0是尺度参数，ξ是形状参数。后验推断遵循贝叶斯定理，即（2）p(θ∣D)∝L(D∣θ)p(θ)，其中θ=(μ,σ,ξ)。为了在样本量有限（N=27）的情况下规范估计，指定了弱信息先验。具体来说，μ～N(x?,2σ^)，logσ～N(0,1)，ξ～N(0,0.5)，其中x?和σ^是年最大降水量序列的样本均值和标准差。ξ的先验以零为中心，具有适度的扩散，这与短时间水文气象记录的常见做法一致[3]，[4]。后验样本使用自适应Metropolis-Hastings算法抽取，共4000个后热抽样（每个链1000个热抽样）。使用split-R^统计量和批量有效样本量（ESS）评估收敛性；所有GEV参数的R^≤1.014且批量ESS≥346，表明混合效果良好。对于返还期T（超越概率p=1/T），计算每个后验抽样θ(m)的返还水平，即（3）zT(m)=F^(-1)(1-p;θ(m))，从而得到设计降雨量的后验分布。

**3.2 基于CRPS的BMA设计降雨量**
使用CRPS根据历史表现确定模型权重，CRPS衡量预测分布与观测值之间的差异，即（4）CRPS(F,x)=∫_(-∞)^∞F(z)-I(z≥x)dz，其中I(?)是指示函数。通过softmax映射获得权重，即（5）wk=exp^(-CRPSk/τ)∑_j=1^Kexp^(-CRPSj/τ)，∑_k=1^Kwk=1，其中k索引模型，基线分析使用τ=20。

**3.3 水文和淹没模型验证**
RRI模型已在[12]中针对卡加延河流域进行了校准和验证，针对的是2020年11月台风Ulysses期间图格加劳（Buntun Bridge）的水位水位曲线。流向数据来自HydroRIVERS[20]，地形数据来自HydroSHEDS；河流几何形状使用卫星图像和高程数据进行了校准。该模型很好地再现了观测到的峰值水位（模拟值为13.28米，实际值为13.30米；NSE=0.72），这对于推导降雨-阶段关系器是足够的。模拟的淹没范围也与2020年11月13日获取的Sentinel-1 SAR图像一致[11]。

**3.4 洪水模拟和降雨-阶段关系**
使用RRI模型进行了洪水模拟，空间时间降雨结构固定为2009年9月至10月观测到的Ondoy–Pepeng事件（7天累积窗口）。该模板在流域平均7天总降水量P=200–800毫米的范围内进行了缩放，并通过模型进行模拟。在图格加劳（Buntun Bridge）分析了水位。拟合了线性对数关系式，即（6）H=aln(P)+b，其中P是流域平均7天降水量（毫米），H是峰值水位（米）。

**3.5 损失估计和降雨-损失函数**
通过将淹没深度与暴露层和深度-损失函数相结合来估计直接损失[11]。对于实验中的每个降雨量P，计算了建筑物、水稻和玉米的事件损失，并在整个研究区域内进行了汇总。对于每种资产类别c，损失通过二次损失函数近似表示，即（7）Dc(P)=α_cP^2+β_cP+γ_c，其中Dc(P)是事件损失（菲律宾比索）。负值被截剪为零。

**3.7 不确定性传播和总结指标**
对于每个情景，从情景条件分布中抽取了100年流域平均7天总降水量的一系列蒙特卡洛样本。每个样本使用公式（6）映射到峰值阶段，使用公式（7）映射到事件损失。使用中心50%和95%可信范围的宽度总结了损失的后验不确定性，即（8）W50(c,s)=Q0.75(Dc,s)-Q0.25(Dc,s)，（9）W95(c,s)=Q0.975(Dc,s)-Q0.025(Dc,s)，其中Q_p(?)表示p分位数，c表示资产类别，s表示情景。此处，W50表示典型结果的分布范围，而W95表示包含尾部的不确定性。

**4. 结果**
**4.1 观测到的极端降雨**
图1显示了流域平均（A）1天和（B）7天年最大降量的返还期图。黑色圆圈表示使用Gringorten公式（10）Fi=i?0.44n+0.12计算出的观测值的实证绘制位置，其中i是按升序排列的排名，n是样本大小。红色曲线显示了来自拟合GEV分布的理论返还水平。拟合的GEV参数分别为1天最大降量的θ={μ=72.4450,σ=22.3369,ξ=?0.0678}，7天最大降量的θ={μ=160.9968,σ=30.9120,ξ=?0.2216}。表1显示，7天的返还水平大约是1天对应值的2.1–2.5倍，且在较长的返还期这个比例更大；100年的7天总降水量达到410毫米，而1天总降水量为161毫米。7天序列的正形状参数表明，与1天序列相比，其尾部不确定性更大，这意味着罕见事件降雨量的估计对多日尺度上的尾部不确定性更为敏感。

**下载：** 下载高分辨率图像（267KB）
**下载：** 下载完整尺寸图像

**图1.** （A）不同返还期的1天降水量总量。（B）不同返还期的7天降水量总量。
**表1.** 不同返还期的河盆平均降水量（毫米）

**4.2 基于BMA的设计降雨量不确定性**
表2总结了基于CRPS的7天年最大降水量的BMA权重。最大的权重分配给EC-Earth3-Veg-LR（w=0.50；CRPS =39.65毫米），其次是CanESM5（w=0.32；CRPS =48.70毫米）和ACCESS-CM2（w=0.18；CRPS =60.79毫米）。拟合的GEV形状参数也因模型而异：EC-Earth3-Veg-LR显示出轻微的重尾行为（ξ?=0.11），而CanESM5和ACCESS-CM2则表明尾部有界限的趋势（ξ?<0）。图2显示了在SSP1–2.6和SSP3–7.0下，基于CRPS加权的BMA后验返回水平曲线与基线的对比。在后两种情景中，后验不确定性随着事件发生频率的降低而扩大，但SSP3–7.0显示出更大的上升趋势和更宽的可信区间。图3展示了RL100的后验分布。SSP1–2.6的分布与基线接近，而SSP3–7.0的分布明显向右移动，尾部更宽。对于分析性传播，RL100降雨样本被近似为GEV分布，其参数分别为：(11)fBMASSP126(z)=fz;μ=407.06,σ=45.50,ξ=?0.0156; (12)fBMASSP370(z)=fz;μ=469.29,σ=69.72,ξ=?0.00386。

表2列出了基于CRPS的BMA权重，用于7天年最大降水量序列。CRPS以毫米（7天总量）为单位报告。

GCM | 权重 wCRPS (mm) | 尾部形状（ξ?的符号）
----|------------|-------------------
EC-Earth3-Veg-LR | 0.50 | ξ?>0 |
CanESM5 | 0.32 | ξ?<0 |
ACCESS-CM2 | 0.18 | ξ?<0 |

下载：高分辨率图像（353KB）
下载：全尺寸图像

图2：7天年最大降水量的GEV概率图，显示了在(A) SSP1–2.6和(B) SSP3–7.0下的BMA后验（红色）。实线表示中位数；阴影带表示95%的可信区间。

下载：高分辨率图像（224KB）
下载：全尺寸图像

图3：在SSP1–2.6和SSP3–7.0下，RL100（100年7天Q100）的BMA后验分布。红色虚线表示基线RL100。

4.3. 对RL100不确定性的洪水阶段响应
图4显示了在SSP1–2.6和SSP3–7.0下，Tuguegarao（Buntun Bridge）的水位水位图。在SSP1–2.6下，后验中位水位图与基线接近，不确定性带较窄；而在SSP3–7.0下，水位图向上移动，不确定性带较宽。

拟合的降雨-峰值阶段模拟器与RRI模拟非常吻合，其关系为：(13)H=6.56ln(R)?27.46，R2=0.99，其中R是流域平均7天降雨总量（毫米），H是Tuguegarao对应的峰值水位（米）。使用这个模拟器，RL100降雨样本被映射到后验峰值阶段分布。图5显示，SSP3–7.0产生的峰值阶段分布比SSP1–2.6的分布向更高值偏移，并且分布范围更广。红色虚线表示在当前气候RL100设计降雨下的基线峰值阶段。

4.4. 损失函数和RL100损失分布
图6显示了从RRI模拟中得出的降雨-损失关系。二次损失函数对所有三类资产都提供了极好的拟合（表3），R2值在0.997到0.999之间。

图7显示了在SSP1–2.6和SSP3–7.0下RL100损失的后验概率密度函数（PDF）。在所有资产类别中，SSP3–7.0将分布向右移动并扩大了分布范围；所有分布都呈右偏态，其中建筑物的绝对损失最大。红色虚线表示使用后验中位数RL100降雨量估算的确定性损失，表明点估计低估了尾部风险。

表3：方程（7）中的二次损失函数系数。Dc(P)以菲律宾比索（PHP）为单位，P是流域平均7天降雨总量（毫米）。

表4：基于RRI模拟，在P=200–800毫米（50毫米间隔）的条件下，建筑、水稻和玉米的流域平均7天降雨总量与事件损失之间的关系，以及[11]中的损失评估方法。实线表示用于不确定性传播的二次拟合。

图7：在SSP1–2.6（左）和SSP3–7.0（右）下，Tuguegarao（Buntun Bridge）的RRI基础水位图。红线：后验中位数RL100响应。阴影带：由第2.5和第97.5百分位数RL100降雨量驱动的水位图之间的范围。黑线：由当前气候RL100设计降雨量驱动的基线水位图。

4.5. RL100损失的不确定性宽度
图8使用W50和W95指标总结了RL100损失的后验不确定性宽度。在所有资产类别中，SSP3–7.0下，中心50%的可信范围宽度W50增加了。建筑物的W50增加了最多，从5.63亿增加到10.76亿PHP，其次是水稻（3.96亿增加到7.2亿PHP）和玉米（3.75亿增加到4.91亿PHP）。在尾部，W95的增加更为明显，建筑物的W95从1.944亿增加到3.817亿PHP。水稻也从13.36亿增加到24.94亿PHP。玉米的W95从9.13亿增加到13.55亿PHP，尽管其较低尾部可能因为负多项式值被截断为零。

5. 讨论
5.1. RL100不确定性的主要影响
当RL100降雨的不确定性传播到洪水阶段和损失时，这种不确定性会非线性放大，其中SSP3–7.0的放大效应明显强于SSP1–2.6。这意味着RL100的点估计可能会低估高影响结果，特别是对于建筑物来说，其绝对不确定性增长最大。因此，对于传达未来的洪水风险，概率报告比确定性RL100估计更为合适。

5.2. 风暴模板的敏感性
为了测试损失估计对降雨结构的敏感性，除了使用Ondoy–Pepeng模板外，还使用了2020年11月的Typhoon Ulysses模板进行了降雨-损失模拟。图9比较了两个模板在每小时流域平均降雨量和7天总量的标准化空间分布方面的差异，而表4和表5总结了由此产生的损失和不确定性宽度。在两个模板下，总事件损失在校准范围内相差约±10%，总损失W95仅相差3%–8%。这些结果表明，主要结论——即在SSP3–7.0下尾部风险的放大——对于这两种代表性事件的选择是稳健的。同时，不同资产类别之间的损失分布更加敏感，这反映了两个模板在降雨时间和空间分布上的差异。

5.3. 基于情景的稳定性和非稳定性考虑
本研究采用基于情景的稳定性方法：分别对历史（1995–2014年）和未来（2040–2099年）模拟期间的年最大序列拟合独立的平稳GEV模型，通过不同时期GEV参数的差异来捕捉气候变化信号，而不是通过显式的时间变化模型。这种设计选择简化了解释，并减少了从可用数据中需要估计的参数数量，但它意味着在2040–2099年期间的趋势和连续的气候演变没有得到解决。

5.4. 稳健性和方法论范围
表6总结了三项稳健性检查，以评估主要结果对关键建模选择的敏感性。所有检查都确认了主要结论——即在SSP3–7.0下尾部风险的放大是稳健的，同时也指出了BMA权重依赖于偏差校正的方法学注意事项。

5.5. 限制和未来工作
- 观测记录较短：观察到的年最大序列仅覆盖了27年（N=27），相对于估计的100年回归周期来说时间较短。贝叶斯推断通过使用信息量较少的先验来减轻但无法消除这种外推不确定性，通过规范尾部估计来缓解这一问题。为RL100报告的宽后验可信区间直接反映了这一观测限制，应将其解释为基于情景的条件概率范围，而不是精确预测。在像CRB这样的数据稀缺的流域中，这种概率方法尤为有价值，因为这些流域缺乏长期的水文记录，明确量化尾部不确定性比隐含它的确定性点估计更具信息性。
- 单一风暴模板：模拟器是通过缩放单一风暴模板（Ondoy–Pepeng）得出的。如第5.2节所示，总体结果是稳健的，但每类资产的损失更依赖于模板（例如，水稻的W95相差高达24%）。测试包括具有不同空间梯度和持续时间比的合成设计风暴的其他模板是未来工作的重点。
- 代理模型校准范围和外推：模拟器在校准范围内基于P=200–800毫米进行校准。根据BMA RL100后验，超过800毫米的样本比例在SSP1–2.6下可以忽略不计（<0.01%），而在SSP3–7.0下大约为0.8%（第99百分位RL100约为787毫米，略低于上限校准值）。因此，对于中心推断来说外推风险较小，但在SSP3–7.0的极端右侧尾部则不可忽略。建议进行额外的RRI模拟，以扩展和验证这些模拟器。
- 损失和脆弱性模型的不确定性：报告的可信范围主要反映了由降雨驱动的不确定性，并没有完全传播深度-损失曲线、暴露值或资产估值的不确定性。此处使用的深度-损失函数来自菲律宾背景下的验证来源[14],[15],[16]；尽管如此，深度-损失参数的±20%变化表明这一成分可能对总损失变化贡献了额外的10%–20%。
- 流域平均降雨量和空间变异性：流域平均与GCM数据的空间分辨率（0.25°）一致，提高了年最大序列的统计稳定性。然而，极端降雨的空间变异性可能会影响淹没模式和损失。Ondoy–Pepeng风暴模板部分地通过保持流域内真实的降雨梯度来补偿这种简化。其次，这种降雨不确定性在灾害阶段会以非线性方式加剧：在SSP1–2.6和SSP3–7.0的情景下，建筑损失95%可信区间宽度（W95）从19.44亿菲律宾比索增加到38.17亿菲律宾比索，几乎翻了一番。第三，使用2020年的台风“尤利西斯”进行的风暴模板敏感性分析表明，综合概率结论对所选风暴模板具有鲁棒性：总损失估计值在±10%范围内一致，而总损失的W95在不同模板之间的差异仅为3%–8%。这些结果具有直接的政策意义。在SSP3–7.0情景下尾部不确定性的加倍意味着，排放路径不仅显著影响預期的洪水损失，还影响成本效益分析的可靠性，因为在更高强迫情况下确定性设计值会低估尾部风险。基于概率可信区间的决策框架（如与仙台框架对齐的风险指标以及在高度不确定性下的稳健决策）比传统的确定性设计标准更适合这种情况。后验阶段分布还允许推导出特定预警级别的超标概率，为疏散规划和防洪屏障投资提供可行的衡量标准。

**作者贡献声明：**
Yuta Kurihara：撰写、审稿与编辑、原始稿撰写、可视化、验证、软件开发、资源管理、方法论制定、调查、资金获取、正式分析、数据整理、概念构思。
Mamoru Miyamoto：监督、方法论指导。