做气候模式后处理的人十有八九都被同一个问题折磨过模式输出的降雨量、气温和站点观测差得太离谱。均值和方差偏掉就算了更头疼的是整个分布形态都对不上——模式里小雨日太多大雨日太少百分位序列歪得没法直接用。你要是直接拿这种数据去跑水文模型、算干旱指数结果基本没法看。Quantile Mapping分位数映射简称QM就是专门干这个事儿的后处理手段。它做的事情本质上是一句话把模拟值的概率分布搬到观测值的概率分布上去。听起来不复杂但实际做起来牵扯到经验分布估计、零值处理、外推控制、训练期与验证期切分等一堆细节。这篇东西我打算把这些年用R语言做分位数映射后处理的实操经验全部倒出来从原理到代码从参数选择到坑点记录尽量写得能直接照着跑。适合谁来读如果你正在做气候模式数据降尺度、气象水文预报订正、或者任何涉及“模式输出需要对齐观测”的活儿这篇文章应该能帮你省下不少试错时间。已经熟悉原理的人可以直接跳到第3节的R语言实现和第5节的问题排查。1. 模式数据差在哪Quantile Mapping补什么1.1 模式系统性偏差的三张脸均值、方差和分布形态先说清楚我们要解决的问题。气候模式从全球气候模式到区域气候模式再到天气预报模式输出的变量通常都带有系统性偏差。这种偏差不是随机噪声而是稳定存在的、有统计规律的误差主要表现在三个层面。第一层是均值偏差。比如某个模式在你研究的区域夏季降水常年比观测偏高30%或者2米气温普遍偏低2℃。这是最直观的偏差很多人一开始只做这层校正。第二层是方差偏差。模式模拟的变量波动幅度经常比实测窄尤其是降水这种高变率变量。日降水量的方差如果被压缩意味着极端事件暴雨、干旱在模式里表现不出来。只调整均值解决不了这个问题。第三层是分布形态偏差。这一层最隐蔽也最麻烦。模式不仅均值、方差不对可能连整体分布的形状都跟观测对不上——偏度不同、峰度不同、甚至右尾拖多长都不一样。比如很多模式在夏季降水上表现为“小雨日太多、暴雨太少”这种偏差用均值校正根本无法消除。1.2 为什么线性缩放救不了分布形态而分位数映射能常见的简单后处理手段比如线性缩放Linear Scaling和方差校正Variance Scaling本质上是做一阶或二阶矩的匹配。线性缩放就是给模拟值乘个系数让平均值对齐方差校正在此基础上再拉伸一下序列让方差也对齐。问题在于这两种方法都预设了“模拟值的分布形态和观测值一致只是位置和尺度不同”。可实际情况中这个假设经常不成立。极端降水、气温极值这些尾部行为偏差模式跟主体区域完全不同——你可能均值只偏了10%但99分位偏了50%。线性缩放对主体区域的校正有效对尾部要么没校到位要么直接过头。分位数映射的思路完全不同它不假设两个分布形态相似而是直接逐分位点建立映射关系把模拟分布的每一个分位点都搬到观测分布对应分位点的位置上。这样校出来的结果不仅均值、方差接近观测整个分布形态——包括偏度、峰度、尾部厚度——都能对齐。这才是QM在气象水文后处理里被当成主流方法的原因。2. 分位数映射的原理一句话版本和数学版本2.1 核心直觉把模拟值的分位点搬到观测分布上去分位数映射的核心思想用一张图就能说清楚把模拟值投影到它自己的累积分布函数CDF上得到对应的概率位置然后拿着这个概率位置去观测分布的CDF反函数上查对应的值。输出值就是校正后的结果。公式长这样sim_corrected F_obs^{-1}(F_sim(sim_raw))其中 F_sim 是模拟值的累积分布函数F_obs^{-1} 是观测值的分位数函数累计分布函数的逆函数。整个过程可以拆成三步算出模拟值在模拟分布里的累计概率比如0.85。拿着0.85这个概率去观测分布里找它对应的分位数值。用观测分布的那个分位数值替换掉原始模拟值。这样做的结果就是校正后的模拟值在观测分布中的相对位置跟原始模拟值在模拟分布中的相对位置保持一致。换句话说模式里最干旱的那一天校正后对应观测里最干旱的那一天模式里百年一遇的暴雨校正后对应观测里百年一遇的暴雨。位置对应关系完整保留这是分位数映射最迷人的地方。2.2 从累计概率到分位数ecdf和quantile这对组合拳在R语言里这个原理落地非常直接。经验累积分布函数用ecdf()一步就能算出来分位数函数用quantile()也能一行搞定。# 构造模拟的历史降水序列假设和观测同期 sim_his - rgamma(1000, shape 2, rate 0.5) # 观测降水序列 obs - rgamma(1000, shape 3, rate 0.4) # 第一步ecdf 得到模拟值的累积概率 F_sim - ecdf(sim_his) u - F_sim(6) # 模拟值6对应的累计概率大概0.8左右 print(u) # 第二步quantile 反查观测分布中的分位数值 corrected_value - quantile(obs, probs u, type 8, names FALSE) print(corrected_value)ecdf()返回的本质上是一个阶梯函数输入一个值输出这个值在样本中的累积概率。quantile()是这个过程的逆操作输入一个概率输出样本分位数。两者配合使用就是最朴素的分位数映射实现。有一个细节值得多说一句quantile()函数的type参数很多人不重视但实际影响很大。R语言里quantile有9种算法默认是type7而水文气象后处理里常用type8因为它在正态分布下是无偏估计。你用的type不同极端分位点的数值可能差不少。一般建议固定用type8保持一致性。2.3 降水的零值堆积必须单独处理降水数据和气温数据有个本质区别降水序列里有大量零值。晴天的降水是0小雨可能也是0.1毫米以下被记成0。这就导致降水的经验分布函数在0那个位置有一个巨大的跳跃。如果你直接对整个降水序列做分位数映射会遇到两个问题。第一大量的零值会让映射曲线在低分位区域变得极度陡峭第二模式模拟的降水日数和观测可能本就不同——模式可能降水日过多你校正完后不仅量不对降水频率也不对。常见的解法是引入“湿日阈值”wet-day threshold。先把低于阈值的日子归为干日高于等于阈值的日子归为湿日然后只对湿日的量值做分位数映射最后把干日全部置0或者按需要保留一些极小值。wet_day_threshold - 0.1 # 单位按你的数据来一般降水用0.1mm或1mm wet_index_obs - obs wet_day_threshold wet_index_sim_his - sim_his wet_day_threshold阈值的选取没有绝对标准一般根据观测数据的精度来定。如果观测降水精确到0.1mm阈值就用0.1mm如果有些自动站把微量降水记成0就用1mm更稳妥。另外这种零值处理方式也意味着QM做出来的校正结果在降水日数上也会向观测靠拢这本身算是一个额外的好处。3. R语言实操从数据准备到映射完成3.1 数据准备三个数据集该按什么规则摆动手写代码之前先把数据结构理清楚。做分位数映射至少需要三个时间序列观测数据obs训练期的站点观测值。比如1981-2010年30年的逐日降水。模式历史模拟sim_his与观测同时期的模式输出。模式未来或预报期数据sim_fut待校正的目标数据。关键约束是观测和模式历史模拟必须在同一时间段上配对因为我们假设两者的经验分布是在同一个气候状态下估计的。而待校正数据可以是任何一个时期只要确保气候背景不出现极端突变就行。实际操作里经常遇到观测数据长度和模式历史模拟长度不一致的情况。比如观测只有25年模式历史有30年。我的习惯是先把两边裁剪到共同的时间段按日期对齐再去训练映射关系。两者不等长的情况下做的ecdf权重会被拉偏映射关系就不准。如果你用的是CMIP6这类全球模式输出建议先插值到站点或目标格点。插值方法对后期QM结果影响不大因为分位数映射处理的是分布而不是空间位置但你要是连格点都没对齐就做映射统计上没问题空间上会被人挑毛病。3.2 核心代码写一个能处理湿日阈值的QM函数直接给出一个我平时常用的经验分位数映射函数这个版本同时处理了降水零值问题也支持对任意长度序列做校正qm_correct - function(obs, sim_his, sim_fut, wet_day 0.1, type 8) { # obs: 训练期观测与sim_his同期 # sim_his: 训练期模式模拟 # sim_fut: 待校正的模式数据任意长度 # wet_day: 湿日阈值气温场景设-Inf即可 # type: quantile函数的分位数算法类型 # 湿日筛选降水场景 wet_obs - obs wet_day wet_his - sim_his wet_day obs_wet - obs[wet_obs] sim_his_wet - sim_his[wet_his] # 湿日的经验分布 F_sim_wet - ecdf(sim_his_wet) F_obs_inv_wet - function(u) quantile(obs_wet, probs u, type type, names FALSE) # 对sim_fut的逐日校正 corrected - rep(0, length(sim_fut)) wet_fut - sim_fut wet_day if (any(wet_fut)) { u - F_sim_wet(sim_fut[wet_fut]) corrected[wet_fut] - F_obs_inv_wet(u) } return(corrected) }这个函数的核心流程是先按阈值分开干湿日对湿日序列做ecdf然后把未来模拟的湿日值映射到观测湿日分位数上干日直接置0。如果你处理的是气温直接把wet_day -Inf传进去这样所有值都会进入湿日序列相当于全程映射。注意气温场景下我一般不建议把干日置0所以要设置一个极大负数阈值确保所有值都保留下来。用小规模数据测试一下函数是否能跑通set.seed(42) obs_daily - rgamma(365 * 10, shape 3, rate 0.5) # 10年观测 sim_his_daily - rgamma(365 * 10, shape 2, rate 0.35) sim_fut_daily - rgamma(365 * 50, shape 2.2, rate 0.4) corrected_fut - qm_correct(obs_daily, sim_his_daily, sim_fut_daily, wet_day 0.1) head(corrected_fut)跑完看一眼结果校正后序列的干日比例、湿日均值应该和观测更接近。这里的gamma分布参数只是示例数据真实项目里直接用观测和模式值替换就行。3.3 用图和统计指标验证校正效果代码跑通只是第一步。验证校正是否有效是后处理里绝不能省的一环。我一般做三件事。第一画分布对比图。把观测、原始模拟、校正后模拟的密度曲线画在一起直接看分布形态是否对齐。这个最直观校正前模拟的分布可能偏左、偏矮、偏胖校正后应该跟观测的密度曲线基本重合。plot(density(obs_daily, na.rm TRUE), col black, lwd 2, main Distribution Comparison, xlab Precipitation) lines(density(sim_his_daily, na.rm TRUE), col red, lwd 2) lines(density(corrected_fut, na.rm TRUE), col blue, lwd 2) legend(topright, legend c(Obs, Sim Raw, Sim Corrected), col c(black, red, blue), lwd 2)第二画QQ图。把观测分位数和校正后模拟分位数做点对点对比如果点基本落在1:1线上说明分布对齐得很好。qqplot(obs_daily, corrected_fut, main QQ Plot: Corrected vs Obs) abline(0, 1, col red, lwd 2)第三算统计指标。我最常用的是RMSE、KGE和KS检验。RMSE看整体误差KGE看相关、均值偏差和变率偏差三个维度KS检验看分布是否统计显著差异。rmse_raw - sqrt(mean((obs_daily - sim_his_daily)^2)) rmse_corr - sqrt(mean((obs_daily - corrected_fut[1:length(obs_daily)])^2)) cat(RMSE raw:, rmse_raw, RMSE corrected:, rmse_corr, \n) ks_raw - ks.test(obs_daily, sim_his_daily) ks_corr - ks.test(obs_daily, corrected_fut[1:length(obs_daily)]) cat(KS p-value raw:, ks_raw$p.value, KS p-value corrected:, ks_corr$p.value, \n)注意一点这里比较校正后模拟和观测必须是同一个时间段。上面示例里corrected_fut是最初假设的未来期长度跟obs_daily不一样所以我用了[1:length(obs_daily)]来截取。真实项目里做验证要用独立验证期的校正结果不能把训练期数据拿来检验这是很多初学的人容易犯的错误。4. 进阶参数映射、DQM和现成R包怎么选4.1 参数化映射降水和气温该配什么分布经验分位数映射的优点是不假设分布形式数据长什么样就映射成什么样。但它的缺点是观测样本尾部的信息有限如果要校正的值超出了训练期样本范围外推情况经验分布就不好使了。这时候参数化分位数映射就有用场。它先假设变量服从某个理论分布然后通过拟合得到分布参数再做映射。参数化映射的优势在于外推区域的表现更平滑、更符合理论预期不会出现“样本外没数据所以映射不出来”的尴尬。降水和气温适配的分布不同。降水是偏态且非负的一般用Gamma分布或者Weibull分布来拟合。Gamma分布有两个参数形状和尺度能较好地描述降水的主体特征。气温比较接近正态分布用正态分布拟合就够了。当然极端高温场景可以考虑GEV广义极值分布但日常使用中正态假设基本够用。R里做参数化QM可以自己用fitdistrplus包拟合分布也可以直接上qmap包后者把这些都封装好了。library(qmap) # 降水用Gamma分布做参数化映射 qm_fit_gamma - fitQmapPAR(obs obs_daily, mod sim_his_daily, qstep 0.01, type P99) sim_corrected_gamma - doQmapPAR(sim_fut_daily, qm_fit_gamma)fitQmapPAR支持多种分布类型降水场景用type P99表示仅校正到99分位超过99分位的值用最大校正因子继续外推。这个细节很关键后面讲外推问题的时候会再展开。4.2 DQM的思想在去趋势之后再做映射传统QM有一个假设训练期建立的映射关系在未来时段依然成立。这在气候平稳时期问题不大但如果你处理的是气候变化情景下的未来数据模式输出的均值可能已经发生了趋势性变化直接用历史和观测建立的关系去校正未来数据会把这个趋势抹掉或者扭曲。Detrended Quantile MappingDQM去趋势分位数映射就是针对这个问题做的改进。它的思想分三步先把未来模拟值减去一个平滑的趋势项得到一个“去趋势”的序列然后用历史时期建立的映射关系对这个去趋势序列做校正最后再把趋势加回去。这样既保持了QM的分布校正能力又保留了模式模拟出来的气候变化信号。换个角度理解QM负责把分布形态校准到观测水平趋势项则代表模式对未来的变化预测。两者各管一摊互不干扰。DQM具体的实现代码稍微长一些核心在于用一个滑动窗口或LOESS平滑来提取趋势再去趋势、映射、加回趋势。MBC包里的mba.correction函数就实现了包括DQM在内的多种去趋势变体如果做气候变化影响评估建议直接研究这个包。4.3 现成工具qmap、MBC、hyfo快速上手自己手写QM函数有个好处是逻辑完全透明出了问题好排查。但在生产环境里我更推荐在成熟包的基础上做二次开发。这里列三个我常用的R包适用场景各不相同。qmap是经典中的经典。fitQmapQUANT做经验分位数映射fitQmapPAR做参数化映射fitQmapSSPLINE做样条平滑映射。doQmap*系列函数统一执行映射。这个包适合大多数常规后处理任务。# qmap包的经验分位数映射对应我自己写的函数 qm_fit_quant - fitQmapQUANT(obs obs_daily, mod sim_his_daily, qstep 0.01, wet.day TRUE, type tricube) sim_corrected_quant - doQmapQUANT(sim_fut_daily, qm_fit_quant)MBC包做的是多变量偏差校正。当你要同时校正多个变量比如降水、气温、辐射并且要保留变量之间的相关性结构时单变量的QM就不够了。MBC包里的函数能做多变量版本的QM复杂度高不少但思路仍然是基于分位数映射的框架。hyfo包主打气象水文领域里面不仅有分位数映射还整合了聚合、插值、可视化等功能适合想在一个包内完成整条后处理流程的人。三个包怎么选我的建议是只想快速校正单变量用qmap做气候变化影响研究且涉及多个变量协变用MBC做水文气象业务预报订正并且想顺便出图用hyfo。5. 后处理路上我踩过的坑5.1 外推区域的极端值爆炸怎么限制用经验分位数映射的时候最棘手的问题之一就是外推。比如观测期30年样本里最大日降水量是200mm。模式未来模拟里的某个格点某一天报了280mm。这个值在模拟分布里可能已经超过了历史模拟范围对应的累积概率接近1但还没到1。这时候用观测的经验分位数函数反查查出来的可能是观测最大值200mm附近的值甚至因为分位数算法在最末端的不稳定直接跳到比200mm还夸张的值。换句话说经验QM天然会把超出训练范围的值往观测极值上压导致尾部的变率被压缩。模式里的极端事件越极端校正后反而越趋同。我常用的解法是给映射加一个限制。一种做法是像fitQmapQUANT里type tricube那样对分位数映射曲线做平滑尾部不再完全按经验分位数走而是按一个局部回归的拓展趋势走另一种做法是设定一个尾部截断点超过某个分位比如99%的值不再继续按分位数映射而是按最大校正因子等比例缩放。# 手动限制外推超过训练范围的值按极值比例缩放 max_obs - max(obs_daily) max_sim_his - max(sim_his_daily) scale_factor - max_obs / max_sim_his # 校正后如果发现值超出观测最大值按比例缩回 corrected_fut[corrected_fut max_obs] - corrected_fut[corrected_fut max_obs] * scale_factor这个做法比较粗糙但胜在简单有效防止校正结果出现明显不合理的物理数值。5.2 湿日阈值选错降水天数直接乱掉湿日阈值看着是个小参数实际上影响很大。阈值设得太低会把大量微量降水0.01mm、0.05mm当成湿日纳入映射而这些微量降水在观测里可能根本没有记录阈值设得太高又会把一部分真实降雨日归为干日导致校正后的降水日数偏少。我自己的习惯是先做数据探索性分析。画一个干湿日比例的敏感性曲线横轴是阈值从0到2mm纵轴是观测和模式模拟的湿日比例差。选一个两边差异最小的阈值或者选一个符合观测数据精度、同时湿日比例差可接受的阈值。比如很多自动站数据精度是0.1mm那就先试0.1mm如果模式的微量降水日太多就逐步提高到0.5mm或者1mm。还有一个小细节湿日阈值的判断训练期和未来期要一致。你按0.1mm训练映射关系未来期也按0.1mm判定干湿日不能换。否则校正后的未来序列干湿日频率会受到两套标准影响结论解释起来很麻烦。5.3 训练期与验证期数据串扰检验结果虚高这是做后处理最容易犯、也最隐蔽的错误拿训练期的观测数据同时用来训练映射关系和验证校正效果。这样得到的RMSE、KGE一定好看但这是“开卷考试”的结果代表不了方法在未来期的真实表现。正确的做法是把数据切段。用training split和validation split。比如30年数据前20年用来训练映射关系后10年作为独立验证看校正效果。如果数据量不够也可以用交叉验证或留一法。在气候变化情景里由于未来数据本身无法和观测对应通常做法是用历史观测和模式历史模拟训练关系应用到模式的未来期数据上作为订正结果然后在历史期内做回代检验即把训练好的关系应用到模式历史模拟上与同期观测比较。回代检验有一点要说清楚如果训练期和验证期是同一段数据回代检验的结果会偏乐观。我一般用交替验证——拿1951-1980年训练1981-2010年验证然后反过来再跑一次。两个方向的验证结果如果都稳定才说明映射关系是可靠的。再补一个容易忽略的点降水序列的时空变率非常大30年训练期可能只覆盖了气候态的几个干湿阶段样本量在尾部仍然不够。这种情况下做经验QM的验证时要特别关注高百分位99%、99.9%的误差不能只看整体均值。5.4 非平稳问题气候变化情景下QM可能会失真最后再提一个部分人可能还没意识到的问题——传统QM在气候变化情景下可能低估极端事件的变化。原因很简单。QM的训练关系建立在历史时期模拟和观测的关系上当你把它应用到未来模式数据时它默认“未来的模型偏差和现在一样”。如果模式的未来降水分布整体抬升了QM会把抬升的均值保留下来因为原始值变大了分位点顺带上移但分位数映射用的是历史时期的分位数对应曲线意味着未来新增的超出历史范围的极值会被压缩到历史极值附近导致未来的极端事件被系统性低估。这个问题没有完美解法但有几个缓解策略可以组合使用。第一用DQM这类去趋势方法在映射之前先把趋势提取出来第二把训练期拉长到尽可能覆盖更多的气候变化范围第三对超出历史范围的值采用参数化映射而非经验映射让尾部的延伸更平滑。6. 写在最后分位数映射的方法论框架并不复杂真正费时间的是细节。我做过不少次类似项目最深刻的体会是方法本身只能保证“统计上对齐”不能保证“物理上合理”。一个后处理结果就算RMSE很低、KS检验不显著仍然需要人去判断它的物理意义是否成立——比如校正后的降水是否还符合区域气候特征趋势信号是否被合理保留空间分布是否连续。这些东西代码不会替你检查。所以我的建议是把分位数映射当成一个工具箱里的标配工具而不是万能钥匙。拿到一个新数据集先做探索性分析理解模式偏差的特征再决定用哪种映射变体最后一定要做独立验证。统计指标和物理判断都过关了这个后处理结果才真的能用。