给洪水预报‘纠偏’：手把手教你用Python实现数值降雨预报的线性缩放(LS)与分位数映射(QM)校正

给洪水预报‘纠偏’：Python实战数值降雨预报的线性缩放与分位数映射校正

暴雨预警发出后，洪水预报模型却给出了与实际情况偏差较大的结果——这是许多水文工程师经常遇到的困境。数值降雨预报数据作为洪水模型的关键输入，其准确性直接影响着预警的可靠性。本文将带你用Python实现两种主流校正方法：线性缩放(LS)和分位数映射(QM)，将原始GRAPES预报数据转化为更可靠的输入源。

1. 环境准备与数据加载

在开始校正前，需要配置合适的Python环境。推荐使用Anaconda创建独立环境：

conda create -n rainfall_correction python=3.9 conda activate rainfall_correction pip install numpy pandas scipy matplotlib xarray netCDF4

假设我们有以下数据文件：

• grapes_forecast.nc：GRAPES模式的原始预报数据（NetCDF格式）
• observed_rain.csv：实测降雨数据（CSV格式）

import xarray as xr import pandas as pd # 加载NetCDF格式的预报数据 forecast = xr.open_dataset('grapes_forecast.nc') print(forecast['precipitation'].attrs) # 查看降水变量属性 # 加载CSV格式的观测数据 obs_df = pd.read_csv('observed_rain.csv', parse_dates=['time'], index_col='time')

常见问题处理：

• 若遇到时间维度不匹配，可使用xarray的reindex方法对齐时间轴
• 单位不一致时（如mm/day与kg/m²/s），需进行单位转换：

  forecast['precipitation'] = forecast['precipitation'] * 86400 # kg/m²/s转mm/day

2. 线性缩放(LS)方法实现

线性缩放基于一个简单假设：预报数据与观测数据之间存在稳定的比例关系。这种方法特别适合系统性偏差明显的场景。

2.1 计算月度校正因子

def calculate_ls_factors(forecast, observed, period='M'): """ 计算各月线性缩放校正因子 :param forecast: 预报数据时间序列 :param observed: 观测数据时间序列 :param period: 分组周期（'M'按月） :return: 各月校正因子Series """ # 合并数据并分组计算月均值 combined = pd.DataFrame({ 'forecast': forecast, 'observed': observed }) monthly_means = combined.groupby(pd.Grouper(freq=period)).mean() # 计算校正因子（观测均值/预报均值） factors = monthly_means['observed'] / monthly_means['forecast'] return factors

注意：率定期建议选择3年以上数据，以覆盖不同气候年型

2.2 应用校正因子

def apply_ls_correction(forecast_series, factors): """ 应用LS校正 :param forecast_series: 待校正预报数据 :param factors: 各月校正因子 :return: 校正后数据 """ corrected = forecast_series.copy() for month, factor in factors.items(): mask = corrected.index.month == month corrected[mask] = corrected[mask] * factor return corrected

效果验证：

# 划分率定期和验证期 train_obs = obs_df['2010':'2018'] train_fcst = forecast.sel(time=slice('2010','2018'))['precipitation'] # 计算校正因子 factors = calculate_ls_factors(train_fcst, train_obs) # 应用校正 test_fcst = forecast.sel(time=slice('2019','2020'))['precipitation'] corrected_ls = apply_ls_correction(test_fcst, factors)

3. 分位数映射(QM)方法实现

当误差分布非线性时，QM方法能更好地保持原始数据的统计特性。我们采用Gamma分布拟合降水数据。

3.1 Gamma分布参数估计

from scipy.stats import gamma def fit_gamma_params(series): """ 拟合Gamma分布参数 :param series: 降水数据序列 :return: (shape, loc, scale) """ # 过滤零降水日 nonzero = series[series > 0] a, loc, scale = gamma.fit(nonzero, floc=0) return a, loc, scale

3.2 实施QM校正

def qm_correction(forecast_series, obs_series): """ 执行分位数映射校正 :param forecast_series: 待校正预报序列 :param obs_series: 观测序列 :return: 校正后序列 """ # 拟合Gamma分布参数 fcst_a, fcst_loc, fcst_scale = fit_gamma_params(forecast_series) obs_a, obs_loc, obs_scale = fit_gamma_params(obs_series) # 校正过程 corrected = forecast_series.copy() for i in range(len(forecast_series)): if forecast_series[i] > 0: # 计算预报值的CDF cdf = gamma.cdf(forecast_series[i], a=fcst_a, scale=fcst_scale) # 用观测分布的反函数获取校正值 corrected[i] = gamma.ppf(cdf, a=obs_a, scale=obs_scale) else: corrected[i] = 0 return corrected

关键参数对比：

4. 效果评估与可视化

校正结果的科学评估需要多角度验证：

import matplotlib.pyplot as plt def plot_comparison(original, corrected, observed, title): fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8)) # 时间序列对比 ax1.plot(observed.index, observed, 'k-', label='Observed') ax1.plot(original.index, original, 'r--', label='Original Forecast') ax1.plot(corrected.index, corrected, 'b-.', label='Corrected') ax1.set_ylabel('Rainfall (mm/day)') ax1.legend() # 累积分布对比 for data, style, label in zip( [observed, original, corrected], ['k-', 'r--', 'b-.'], ['Observed', 'Original', 'Corrected']): ax2.plot(np.sort(data), np.linspace(0,1,len(data)), style, label=label) ax2.set_ylabel('CDF') ax2.legend() plt.suptitle(title) plt.tight_layout() return fig

典型评估指标计算：

from sklearn.metrics import mean_squared_error, r2_score def evaluate_performance(original, corrected, observed): metrics = { 'RMSE_original': np.sqrt(mean_squared_error(observed, original)), 'RMSE_corrected': np.sqrt(mean_squared_error(observed, corrected)), 'R2_original': r2_score(observed, original), 'R2_corrected': r2_score(observed, corrected), 'Bias_original': np.mean(original - observed), 'Bias_corrected': np.mean(corrected - observed) } return pd.DataFrame(metrics, index=['Value'])

5. 实战技巧与进阶优化

在实际业务系统中，还需要考虑以下增强措施：

5.1 动态窗口校正

def dynamic_window_correction(data, window_size=365, method='LS'): """ 滑动窗口动态校正 :param data: 包含forecast和observed列的DataFrame :param window_size: 滑动窗口大小（天） :param method: 校正方法（'LS'或'QM'） :return: 校正后序列 """ corrected = [] for i in range(len(data)): start = max(0, i - window_size) window = data.iloc[start:i] if method == 'LS': factor = window['observed'].mean() / window['forecast'].mean() corrected.append(data.iloc[i]['forecast'] * factor) elif method == 'QM': # 简化的滑动QM实现 fcst_params = fit_gamma_params(window['forecast']) obs_params = fit_gamma_params(window['observed']) cdf = gamma.cdf(data.iloc[i]['forecast'], *fcst_params) corrected.append(gamma.ppf(cdf, *obs_params)) return pd.Series(corrected, index=data.index)

5.2 混合校正策略

结合LS和QM的优势：

1. 先用LS校正系统性偏差
2. 对残差应用QM校正局部误差

def hybrid_correction(forecast, observed): # 第一步：LS校正 factors = calculate_ls_factors(forecast, observed) ls_corrected = apply_ls_correction(forecast, factors) # 第二步：QM校正残差 residual = observed - ls_corrected qm_adjustment = qm_correction(forecast - ls_corrected, residual) return ls_corrected + qm_adjustment

性能对比：