作为常年折腾农业水文模型的人,我对SWAP(Soil-Water-Atmosphere-Plant)的感情很复杂。这模型功能确实强,从土壤水运动到作物生长、盐分运移都能模拟,但用起来也真考验人——数据怎么整、参数怎么调、结果怎么分析,每一步都有坑。最近完整跑完一轮“数据制备→敏感性分析→气候变化影响评估”的流程,踩了不少雷也总结了不少经验,整理出来给需要的朋友参考。这套流程无论是写论文、做项目还是评估区域水资源安全,都能直接用上。
1. SWAP模型运行的底层逻辑与数据体系拆解
1.1 模型核心模块与水分运移计算逻辑
想顺利跑通SWAP,先得把它当成一个“土壤水分账户系统”来理解。核心是基于Richards方程的垂向土壤水动力学,一层土一层土地追踪水分进出。模型把作物根系层划分为多个土壤层,每层有自己的水力特性参数,再结合上边界(降水/灌溉/蒸发)和下边界(地下水位/排水),逐日计算含水量和压力水头。
实际计算时,Richards方程本身靠数值迭代求解,这意味着初始条件和边界条件一旦给错,结果直接崩。最常见的是初始水分剖面和实际田间情况差异太大,模型前几十天都在“自我修正”,含水量曲线会跑出明显锯齿。所以我看很多人直接填一个均一含水量做初始条件,这种做法要小心,最好给两层以上不同含水量,靠近地表的稍微干一点,深层接近田间持水量。
另外,模型内部还有作物生长模块和盐分运移模块。作物模块用简单的积温驱动叶面积指数和根深增长,再通过根系吸水函数(Feddes模型)耦合水分胁迫;盐分模块则跟随水流方程求解对流-弥散方程。这些模块单独跑都没事,一耦合就容易出数值震荡,时间步长设置太大时尤为明显。建议先跑通水分模块,稳定了再逐步把作物和盐分模块加进去,这个顺序不能乱。
1.2 模型输入数据的四大类体系全景
SWAP的输入说多不多说少不少,归纳起来就四类:气象数据、土壤水力参数、作物参数、管理措施与边界条件。很多人一开始就被文件格式劝退,其实理解了之后会发现它只是在用文本文件做“数据库”。
气象数据是驱动模型的最核心输入,包含逐日降水量、最高最低温度、太阳辐射或日照时数、水汽压、风速。缺一不可,尤其是辐射和风速,直接决定蒸散量计算的准确性。如果实测站只有降水温度,那就得用FAO推荐的公式估算辐射和ET0,精度能接受,但要注意和站点实测数据之间的系统偏差。
土壤水力参数走的是van Genuchten曲线,包括θr、θs、α、n、Ks这几个。不做实验直接拿文献值也行,但同一套参数换一个地区就可能水土不服。最稳妥的做法是用土壤传递函数(PTF)初步估一套,再用实测含水量剖面微调α和n,这两个参数对结果影响最大。
作物参数最容易忽略但最影响情景对比,比如根深、叶面积指数曲线、作物系数Kc、最大根系吸水深度和胁迫阈值。这些参数是品种和当地农业措施的函数,网上数据库的默认值只能当起点。我的做法是拿当地灌溉试验站的实测蒸发蒸腾量反推Kc曲线,实测拟合下来比任何默认值都靠谱。
管理措施和边界条件主要是灌溉计划、排水基准和水位边界。下边界是模型里最容易出错的地方,自由排水边界、水头边界、还是地下水位边界,选错的后果是底层含水量完全失真。排水设计不好还会导致水分在底部滞留,模拟出来的径流过程线能给你整条河都修出来。
2. 气象与土壤核心数据的制备流程与实操要点
2.1 气象数据格式转换与缺测处理全流程
SWAP的气象文件格式历来是新手第一道坎,它要求一个纯文本文件,数据顺序不能错,而且对日期格式极其敏感。我拿到的原始数据通常是Excel或CSV,第一件工作就是把日期、降水、Tmax、Tmin、辐射、风速、水汽压逐列对齐,写成标准顺序的txt。
数据格式的具体要求是每行一天,字段顺序依次是年月日(一般是4位年+2位月+2位日)、降水量(毫米)、最高温度、最低温度、水汽压(kPa)、风速(m/s)、太阳辐射(MJ/m2)。注意温度和辐射的单位,很多数据源给的是度和小数,不用换算;辐射如果给的是W/m2日均值,要先乘以86400再除以1e6,换算成MJ/m2,这个换算扣掉不是小数——有朋友直接拿W/m2塞进去,算出来的蒸散量偏差超过30%。
缺失数据是绕不开的现实问题,气象站全年总有几天数据缺失。缺测处理我按优先级来:缺口在5天内用相邻日期的线性插值;连续缺一周以上就建立该站的月平均日变化曲线,用多年同期平均值填充;如果整个月都缺,死活都要去附近气象站借数据出来回归,千万别拍脑袋填。插值完一定要画时序图,粗看一眼有没有出现负降水、Tmax小于Tmin这类低级错误。
细心的朋友可能还会遇到降水数据从雨量计换成了自动站导致的不连续问题。这种站点变更引入的系统偏差,常规插值补不了。我的办法是对重叠观测期做线性校正,把自动站数据回归到手动站水平,然后再接入模型。不处理这层偏差的话,年降水量的突变会被模型误判成气候变化信号,非常误导。
2.2 土壤水分特征参数估计与率定实验
土壤参数这块,条件允许的优先用环刀法实测水分特征曲线,然后拟合van Genuchten参数。取土按剖面分层来,表层15cm以内一层,30cm一层,60cm一层,90cm以下可以两层并一层,因为深层土壤的变异主要受母质控制,层位太多反而没有意义。
没有条件做实验的,用PTF函数估算是性价比较高的路径。目前用得比较广的是基于砂粒、粘粒、有机质含量估算θr、θs、α和n的回归公式,比如Rawls等人的PTF模型。但我必须提醒的是,PTF估算值只能做初值用,务必用田间持水量或凋萎含水量的实测点做校正。具体做法是固定θr和θs在合理范围,只调α和n,让模拟的田间持水量匹配实测值。实测和田持之间误差可以控制在5%以内的目标去调,不要追求全曲线完美拟合,因为α和n的敏感性差异很大,一起调容易过拟合出荒谬组合。
饱和导水率Ks是另一个大坑。这个参数在模型里的变幅可以跨好几个数量级,对深层排水和地表径流的模拟结果影响非常大。现场用双环入渗仪测最好,没有条件就用经验公式估算。需要注意的是,即使做实验,尺度效应也很明显——环刀取的小土柱测出来的Ks往往比田间实际值低一两个数量级,因为田间有裂隙、根孔、虫洞这些大孔隙优先流路径。我一般的做法是把实验室Ks乘一个放大系数,具体放大几倍要看当地土壤结构和根系发育情况,没有硬标准,多跑几次对比一下实测排水量就能找到合理尺度。
深层排水和地下水交换这块,很多人不重视,但SWAP最重要的优势之一就是能耦合地下水。如果模拟点位有地下水位观测序列,把下边界设成水位边界,效果比自由排水要好得多。水位边界的输入需要水位动态时间序列,插值到逐日,然后模型会自动解算底层饱和区的水分交换。
3. 敏感性分析的核心方法与系统化实施步骤
3.1 敏感性分析的参数初筛与定性判断
很多人一上来就问用哪个敏感性分析方法,但忘了前面还有关键一步:参数范围与分布怎么定。这个环节对结果影响甚至超过方法本身。范围定窄了,真敏感参数可能被低估;范围定宽了,伪敏感参数又可能被放大。我的经验是每个参数的范围参考两个来源:一个是从文献里收集不同土壤类型或品种的参数上下限,另一个是实测数据的不确定性区间,两个取并集再打九折,作为初步范围。
参数筛选我首推Morris方法,也叫Morris screening或EE法(Elementary Effects)。这个方法需要采样多次,每次只改变一个参数的值,沿着预先设计的轨迹走一遍,最后算出每个参数的均值μ和标准差σ。均值反映参数对输出变量的主效应,标准差反映交互效应。用这个方法跑一遍,可以把十几二十个参数快速分成三类:必须精确定标的、可以先固定的、完全无所谓的。
我自己实际跑过的案例里,土壤参数中α和n基本是最敏感的前两名,Ks和饱和含水量θs紧随其后;作物参数里最大根深和Kc曲线峰值最敏感,而播期之类的管理参数在气象驱动主导的情景下反而不那么敏感。当然这个排序换了地点换了输出变量会变,所以每到一个新站点重新筛一遍才是正途,直接拿文献里的结论套用是不靠谱的。
3.2 定量方差分解与主效应交互效应计算
初筛之后,对筛出来的参数子集再做定量分析,目前主流做法是Sobol方差分解。它能给出每个参数的一阶敏感性指数(主效应)和总效应指数(主效应+交互效应)。
做法简述如下:先生成两批独立的参数样本,第一批作为基准,第二批随机重排。评价指标可以选生育期总蒸散、深层渗漏量或产量模拟值,然后按Sobol公式估算各阶指数。为了保证结果稳定,总计算次数通常需要几百次。SWAP单次运行时间如果小于几秒,完全可以在一台普通电脑上跑完几千次模拟。这比很多区域气候模型的敏感性分析要友好得多,值得利用这个优势多跑几遍,而不是省计算量草草了事。
Sobol分析的代码并不复杂,Python写一套也不难。不过用现成的SALib库确实更快更稳,采样、分析都封装好了,直接调用就行。注意别自己写采样,因为Sobol序列要保证低差异,随机抽样的收敛速度太慢,误差大。
4. 气候变化影响评估的情景设定与批量模拟方案
4.1 未来气候情景数据选择与偏差校正
气候变化影响评估的核心不是跑模型,而是怎么把GCM(全球气候模式)的输出转成SWAP能用的逐日气象数据。这个环节对最终结论的可靠性影响巨大,但它在项目中往往被压缩到最低优先级,实属本末倒置。
GCM数据的获取路径已经很成熟,CMIP6的多个模式都提供各个情景的逐日输出。不过空间分辨率一般都在100公里量级,直接拿来驱动一个点位尺度的SWAP模型,完全不行。这时候需要降尺度。动力降尺度要用区域气候模式,计算成本高;统计降尺度相对轻量,常见的方法包括Delta法、分位数映射法和天气发生器法。
Delta法是简单但容易操作的:用未来时期GCM模拟的月平均变化量(降水比值、温度差值)叠加到历史观测序列上,保持观测序列的波动特征不变。这样做的好处是至少保留了站点尺度的日变异性和极端天气结构,不会因为降尺度过程把日降水分布搞得过于平滑。缺点是假设了变化量在月内均匀,不适合强非线性的变量,比如极端降水。分位数映射法能克服这个缺点,逐月建立历史模拟和观测的分位数映射函数,再应用到未来模拟上,结果一般比Delta法好,但计算量也大一些。
4.2 SSP情景选择与批量模拟的工程化落地
情景选择上,CMIP6的SSP1-2.6、SSP2-4.5、SSP5-8.5基本覆盖了从温和到激进的所有主流路径。农业影响评估我个人推荐多跑几个情景做集合,最少也跑SSP2-4.5和SSP5-8.5两个,前者是“中等发展路径”,后者是“高排放高挑战”。只跑单一情景,审稿人或项目委托方很容易质疑代表性。
批量模拟是另一个新手容易手忙脚乱的环节。SWAP每个情景每套参数都要跑那么多年逐日模拟,手工改输入文件显然不行。我的方案是分三步走:
- 第一步,写一个Python脚本生成所有输入文件组合,把气象数据替换、参数微调、输出路径全部参数化。
- 第二步,调用SWAP的可执行文件逐组运行,每组模拟结束检查日志文件里有没有异常退出或警告信息。
- 第三步,汇总所有输出文件,统一解析成DataFrame格式,按情景、参数组、时间维组织好,后续出图出表时就非常方便。
跑批量之前还有个容易被忽略的工作:确定模拟的预热期。SWAP模拟土壤水分从一个初始状态达到平衡需要时间,一般裸土要跑1到2年才能让水分剖面稳定下来,有作物覆盖的时间更长。所以正式情景分析的时候,建议先跑一段预热期,预热期结果不纳入统计。我有一次忘了设置预热期,直接用初始条件跑未来80年,前10年的土壤水分离群得离谱,整个趋势分析都被带偏,后来花了一周重新筛选处理,教训很深刻。
4.3 气候变化影响指标提取与可视化要点
模型跑完只是第一步,真正有意义的是把原始输出提炼成决策者关心的指标。常用的输出变量包括生育期蒸散量、深层渗漏量、灌溉需水量、相对蒸散胁迫天数、产量或生物量模拟值。把这些变量按年整合成均值、百分位、趋势斜率,再对比基准期与未来时期的变化量或变化率,就完成了气候变化影响的定量表达。
可视化方面建议画三类关键图:一类是年际变化时序加滑动平均,直观展示趋势方向和波动;第二类是空间分布图,如果模拟了多个站点,把变化量填到站点位置做插值;第三类是概率密度转移图,展示未来时期比基准期分布整体偏移的情况,这类图对非专业受众的冲击力很强。
这里重点提醒一句:变化趋势有没有统计显著性一定要检验。很多项目只把均值一减就下结论,但样本量小的情况下均值差异可能根本不显著。做一个简单的Mann-Kendall趋势检验或者两时期T检验,几行代码的事,却能给结论增加不少说服力。
还有一个经验是,气候变化影响评估结果不要只报一个模式的值。GCM本身的不确定性非常大,单模式结果只能算作一个可能的未来。理想的做法是多个GCM做集合,画出集合区间条带图。如果计算条件有限,至少也要选两三个模式取平均,否则结论很难站稳。
5. 全流程常见问题与排查技巧实录
5.1 数据格式、参数超界与迭代不收敛的典型报错
跑SWAP这几年,我积累了一份“报错排查手册”,在这里分享最经典的五类问题以及对应的解决路径。
第一类是气象文件读取出错。症状是模型一启动就报错退出,或者读入的降水量全是零。原因五花八门:文件编码不是ASCII、日期格式不连续、字段之间用了逗号分隔而不是空格、末尾多了一个空行。解决方案是写一个校验脚本,一行一行检查字段数和数值类型,发现异常就定位到具体行,不让模型去“啃”脏数据。
第二类是土壤参数越界。van Genuchten参数不是随便填的,α和n的组合必须保证水分特征曲线单调。如果填的参数产生不了合理曲线,模型在迭代土壤含水量的时候很容易发散。出现这种现象时先检查α是否落在0.001到0.1范围内,n是否在1.05到3之间,θs减θr是否小于孔隙度。超出这个范围就要质疑参数来源了。
第三类是迭代不收敛。SWAP在土壤含水量接近饱和或极干时,Richards方程求解的收敛性会变差。典型表现是模拟到某一天时报错,或者那一时段含水量出现跳变。处理方法有三板斧:降低最大时间步长、减小水含量判断容差、检查上下边界条件是否有突变。从经验看,七成以上的不收敛问题出在边界条件突变,比如突然来了一场大暴雨,而最大时间步长还设置在一天量级,模型根本追不上这个变化。这种情况下,把时间步长在降水日自动加密,问题基本都能解决。
第四类是作物生长异常。如果作物模块报出的LAI或者根深出现负值或者停滞不前,先检查积温计算是否正确。SWAP的积温对不同作物有不同上下限温度,填错温阈会导致积温不增长,作物一辈子长不大。还有一种情况是根深超过了土壤剖面深度,模型会把根延伸到底部边界外面,也会引发异常,需要检查作物参数设置和剖面深度是否匹配。
第五类是水分不平衡。跑完模拟后检查全周期水量平衡,如果总降水量减总蒸散减深层渗漏减土壤储水变化不等于零,模型就有漏洞。这个差值一般来自数值误差,控制在3%以内是正常的,超过5%就要检查是不是有侧向流或回归水没有考虑,或者排水模块参数设置错了。
5.2 不同应用场景下的配置调试策略
针对不同应用目标,SWAP的配置调试策略差异很大,我总结了三类常见场景的配置重点,新手上路时可以直接对号入座。
第一个场景是模拟灌溉农田的耗水规律。这个场景下灌溉策略的设定基本主导了模拟结果。需要重点处理的是灌溉日期和灌溉量的确定,建议在管理措施文件中设置阈值灌溉制度,让模型根据土壤水分亏缺自动触发布水,这样就能分析不同灌溉阈值下的耗水规律。如果设置成固定日期固定水量,模拟出来的结果就只是一个特殊管理方案的表现,推广性有限。
第二个场景是模拟盐分运移和盐碱化趋势。这个场景对初始盐分剖面和底部边界溶质浓度高度敏感。初始盐分建议用实测剖面数据插值,底部边界的盐分浓度不能设为常数,最好给一个季节性波动序列。另外,作物对盐分胁迫的响应参数(阈值EC和斜率)非常关键,直接决定盐分累积到多少开始减产,这个参数在不同作物间差异很大,别拿去其他作物的默认值硬套。
第三个场景是区域尺度的多站点模拟。这种场景最大的问题是参数空间异质性。土壤参数在不同站点间变化很大,如果统一用一套参数跑所有站点,那结果只能反映“平均土壤”的情况,不能反映空间差异。有条件的按土种类型分组率定,每组给一组参数;没条件的也要至少按砂土、壤土、黏土分三类设置,不能一把尺子量到底。
还有一个比较容易被忽略的是时间尺度问题。如果模拟长期变化(比如50年或80年),作物参数和土壤性质需要做动态调整吗?严格来说,土壤有机碳变化可能影响土壤水力性质,但这个模型本身不舍模块去模拟这种长期演变,实际操作中保持参数不变是主流做法,但要在结论中注明这个假设。
5.3 模拟结果合理性检验清单与快速可视化诊断
所有模拟跑完之后,必须做结果合理性诊断再往下分析。我的习惯是画三张诊断图:第一张是土壤含水量时程图,叠加降水事件,检查降水后有没有含水量攀升、干旱期有没有持续下降;第二张是蒸散量对比图,把模拟ET和实测ET0放在一起看趋势;第三张是深层渗漏累积曲线,看年度总量是否和降水变化对应。
如果含水量时程图在降水后没有响应,大概率是入渗参数设置过小,或者表层厚度设置太厚;如果蒸散量常年贴着ET0走,说明作物胁迫没有触发,Kc阈值或根深分布需要调整;如果深层渗漏年际变化和降水脱节,可能底部边界设置的问题太大。
还有一个小技巧是拿实测产量或生物量数据做历史验证。模型模拟的历史产量和统计年鉴的产量做相关分析,R方大于0.6基本可以认为模型有预测能力,小于0.3就要回头找问题了。没有实测产量对比的情况下,至少也要检验生育期内蒸散量的多年均值是否在合理范围——比如华北地区冬小麦全生育期蒸散量一般在450毫米上下浮动,算出个700毫米来肯定有问题。
6. 实战心得与进阶扩展建议
跑完整个流程之后,我最大的体会是SWAP这套东西的瓶颈不在模型本身,而在数据质量和工作流管理。一个模型能不能出靠谱的结果,七八成取决于输入数据的可靠性,剩下三成才看参数率和情景设计。这不是SWAP独有的问题,但在这个模型上体现得特别明显——它对土壤水分特征参数极其敏感,你把参数调好了,复杂情景跑起来水到渠成;反过来,参数不准,后面做了再花哨的敏感性分析和气候变化评估也是空中楼阁。
数据处理和工作流管理方面,几条经验很值得分享。一是所有原始数据不要手动改,全部用脚本处理,保留处理前的原始文件和处理后的标准文件,这样出了任何问题能追溯到底是哪一步引入的。二是文件名里要包含版本号和日期,模型参数文件每次调整都另存一个新版本,别老在一个文件上覆盖,不然跑完一圈下来都不知道结果是哪套参数的产物。三是养成检查日志文件日志文件的习惯,SWAP每轮运行都留有日志,别看它不起眼,里面全是模型运行情况的说明,提取关键警告信息汇总成表,比盯着ctl文件输出轻松得多。
关于扩展方向,除了SWAP自身,有几个组合玩法值得研究。模型和优化算法结合做自动参数率定,目前用PEST或OSPOT这类参数估计工具能自动搜索最优参数组合,省去手动调整的无尽循环。再有就是模型与GIS结合,把单点模拟扩展到田块或小流域,有朋友做出来效果很不错。还可以把SWAP的输出作为作物模型的外部驱动数据,搭一个水分-作物耦合模拟链,这在粮食安全评估项目中很受欢迎。
我现在的工作流已经基本稳定在四个脚本构建的框架内:一个预处理数据、一个批量跑模型、一个做敏感性分析、一个做结果可视化和指标提取。整个流程跑一遍,从拿到站点的原始气象土壤数据到产出全套气候变化影响图表,大约需要一周左右的时间,其中大部分时间花在数据清洗和参数率定上。对刚开始接触SWAP的朋友,我建议不要一上来就追全流程,先拿一个站点的数据,把模型跑通、把输出看懂,再逐步往里面加模块、加情景。模型这个东西,通了流程才会有底气去调参数、去解释结果、去应对审稿人的尖锐问题。祝各位早日跑通自己的第一条SWAP流水线,少走我走过的弯路。