华北平原农田生态系统的Noah-MP模型参数敏感性分析

王 茜，查元源

（武汉大学水资源工程与调度全国重点实验室，武汉 430072）

0 引 言

农田生态系统与人类生命、生活、生产息息相关，它为我们提供赖以生存的粮食，同时也消耗着全球淡水抽取量的70%和总耗水量的90%[1]。作为水-粮食-能源三者相互联系的复杂系统的缩影，准确刻画农田生态系统中涉及的水、能量及CO2通量变化过程，为区域水资源配置、粮食生产、生态环境保护提供理论依据，对保障国家水安全、粮食安全、生态安全有着重要意义。陆面模式作为模拟陆面上发生的各种物理、化学、生物和水文过程，以及这些过程与大气的相互作用的有效工具，历经简单模式阶段与生物大气模式阶段，进入了新一代多模式阶段[2]。从简单地考虑水量平衡，到将陆气间水分与热量的交换过程进行耦合；从不考虑植被的作用，到显式引入植被的生物物理作用，再到通过耦合叶片的光合作用与蒸腾作用实现陆气交换过程的水碳耦合，模式的物理机制逐渐完善，对陆面过程与陆气交换的刻画也逐渐细化。新一代陆面模式中，Noah-MP[3]模型集合多种参数化方案选项，提供多达4 608 种参数化方案组合，方便使用者依据适用性灵活选择，具有广阔的应用前景。

随着陆面模式的不断完善，参数不确定性成为了限制模型模拟性能的主要因素。进入大数据时代，卫星、无人机等多种遥感技术快速发展，使得基于天-空-地多源观测数据建立高分辨率的土壤属性与植被属性等地表参数库成为可能[4]。因此，对陆面模式进行参数敏感性分析，识别出影响模型模拟效果的敏感参数，不仅能够指导参数优化以提升模型模拟性能，也能为地表参数库的构建指引方向。

近年来，许多学者展开了Noah-MP 模型的参数敏感性分析研究。在密西西比河流域评估Noah-MP 陆面模式模拟水文变量的表现时，CAI 等[5]简单地对3 个径流模拟关键参数进行了敏感性测试，并在此基础上以纳什效率系数（NSE）为性能指标进行参数优化，从而提升模型径流模拟能力。ARSENAULT 等[6]利用10 个国际FluxNet 站点的数据对动态植被参数、静态植被参数及土壤参数进行了全局敏感性分析，找出了12 个敏感参数，并指出针对特定区域的具体问题需要对参数分布范围进行进一步约束。而HUO等[7]重点关注草及常绿阔叶林2种典型植被，提出了一种空间抽样与参数抽样相结合的新方法，在美国大陆的网格单元中筛选出植被总初级生产力（GPP）和潜热通量（LE）模拟的敏感参数。SHU 等[8]在黄土高原的研究显示，相比于直接同化叶面积指数（LAI），基于参数敏感性分析的参数优化能更大幅度地降低叶面积指数的模拟值与观测值的均方根误差（RMSE）。上述研究均表明，受不同的下垫面及气候区的影响，需开展针对性的Noah-MP 模型参数敏感性分析试验，以助挖掘参数优化对模型模拟性能提升的巨大潜力。目前，已有大量有关Noah-MP 模型在华北平原农田生态系统的应用型研究[9-11]，但在该区域开展的参数敏感性分析研究未见报道。华北平原是确保全国粮食安全的重要基地[12]，基于Noah-MP 模型开展参数敏感性分析，有助于准确量化这一典型农业生态系统模式与大气间的辐射传输、物质和能量交换，对推动区域经济、社会和生态可持续性协调发展有着促进作用。

1 数据与方法

1.1 观测站点与数据

中国国家生态系统观测研究网络（CERN）的禹城农业综合试验站（116.570 2°E，36.829 0°N，平均海拔22.0 m）与栾城农业生态系统试验站（114.699 3°E，37.893 0°N，平均海拔50.1 m）均位于华北平原，属暖温带半湿润季风气候区，以冬小麦-夏玉米为主要作物，一年2 季。考虑到华北平原轮作特点，即每年10月至次年6月为冬小麦生长季，6月至10月为夏玉米生长季，结合数据质量与可获取性等原因，本研究以2020年6月10日至2021年6月9日的轮作年为研究时段。获取的观测数据包括基本气象数据、多层土壤水（0.1、0.2、0.3、0.4和0.5 m处）及水热碳通量。研究时段内的数据存在大量缺失，需要进行插补处理。插补过程参考R 语言的fluxnetLSM[13]程序包，该程序包主要用于插补遵循Fluxnet 2015 协议的数据。处理后的气象数据仍存在大量连续缺省值难以作为模型所需输入数据，故不加以利用，而土壤水与通量观测数据缺失率低于15%，可以用于模型敏感性分析检验。

1.2 Noah-MP模型

Noah-MP 模型（The Community Noah Land Surface Model with Multi-Parameterization Options）是在Noah 3.0 版本的基础上进行改进和完善的新一代陆面模式。Noah 模型将陆面简单视为0～0.1、0.1～0.4、0.4～1.0 和1.0～2.0 m 的4 层土壤，能量平衡计算在土壤与大气的交界面展开，Noah-MP 模型在此基础上调整模型架构，将冠层与地表分开，引入3 层雪层，同时对积雪/融雪、产汇流、植被生长等关键过程进行优化，并为表1 所示的12 个陆面过程提供多种参数化方案选项。Noah-MP 模型能很好地模拟地表能量收支平衡及陆地水文循环过程，已被作为陆面参数化方案嵌入GFS、CFS 及WRF 等气候/天气预测模式中。为简化输入输出文件读取过程，修改网格化模型的输入输出端口，得到用于站点敏感性分析的单点模型（https://github.com/wangxi1998/NoahMP-Sensitivity-NCP/tree/main/POINT_NOAHMP）。

表1 模型参数化方案设置Tab.1 Model parameterization scheme

1.3 Noah-MP模型参数敏感性分析设计

1.3.1 模型模拟试验的设计

模型所需的气象数据来自GLDAS-2.1[14,15]，空间分辨率为0.25°，时间分辨率为3 h。该数据为Noah-MP 等陆面模式的“在线”模拟集合，故常用作Noah-MP 模型的气象输入数据。利用GEE 平台提取2 个站点在研究时段的数据（2 m 气温、2 m 比湿、10 m 风速、气压、下行短波辐射、下行长波辐射及降水），将其插值得到0.5 h分辨率的数据作为输入。模型初始场参考模型指导手册确定，分层土壤温度与地表温度取全球均值（274.0、278.0、280.0、284.0 和287.0 K），土壤含水率为饱和含水率的40%（取为0.24 m3/m3）。

选用的参数化方案及参数参考ARSENAULT 等[6]和HUO等[7]的研究（见表1 和表2）。选取的42 个参数包含30 个植被相关参数与12 个土壤相关参数；植被相关参数中有12 个与辐射传输过程密切相关，其余参数为动态植被模拟模块涉及参数。模型模拟过程中，模型将土地利用类型与土壤类型作为索引条件，再在模型自带的查找表中匹配对应参数。因此，大部分参数的取值上下限通过查找表确定。少量参数间具有一定大小关系，如冠层顶部高度HVT 与冠层底部高度HVB、田间持水量SMCREF 与孔隙度SMCMAX、凋萎系数SMCWLT与田间持水量SMCREF 及蒸发停止的土壤含水率SMCDRY，需要对它们的取值范围进行约束。

表2 敏感性分析选用的参数及其取值范围Tab.2 Selected parameters and the value ranges for sensitivity analysis

模型运行步长为0.5 h，并利用研究时段内的气象数据进行多轮预热。为了对模型参数敏感性进行系统性的分析，本研究选取潜热通量（LE）、显热通量（H）、净生态系统碳交换量（NEE）、第1 层土壤含水率（SM1）及第2 层土壤含水率（SM2）进行分析。ARSENAULT 等[6]通过数值试验发现当样本量接近1 500 时，全局敏感度基本收敛。更重要的是，那些在小样本量下表现出较高灵敏度的参数在较大样本量下也表现出较高灵敏度，反之亦然。为了减少计算成本，先利用ARSENAULT 等[6]研究中的算例对42 个参数抽样分别抽样250次、500 次、1 000 次与1 500 次，探究合理的抽样次数。有研究表明，动态植被模块开启时，气孔阻抗的土壤水分限制因子的参数化方案对模型参数的取值有较大影响。因此，除对2个站点进行参数敏感性分析外，以CLM 和Noah 2 种气孔阻抗的土壤水分限制因子方案为代表，探索基于土壤水势与基于土壤含水率的参数化方案对参数敏感性的影响。

1.3.2 参数敏感性分析评估方法

本研究采用Matlab 工具箱中的拉丁超立方抽样法（LHS）对参数抽样，抽样得到的参数代入模型运行后，基于Sobol 方差分析计算参数的敏感性指标（实现代码获取网址为https://github.com/wangxi1998/NoahMP-Sensitivity-NCP/tree/main）。全局敏感度Ti（Ti越大代表该模型输出变量对第i个参数越敏感）的计算方法如下[6,16]：

式中：x 为进行敏感性分析的N 个参数的向量（N=42）；xi为第i个参数，x～i为除xi外的N-1个其他参数的向量。蒙特卡洛抽样的集合平均计算方法如下：

式中：M 为拉丁超立方抽样次数；xm为第m 次抽样过程得到的参数向量；是拉丁超立方抽样的2个平行抽样与的参数；f(x)为性能指标。

性能指标选用衡量预测值与真实值之间偏差的均方根误差，计算方法如下：

式中：K 为模型输出变量的时间维度；Rk为第k 个模型运行步长时刻的输出结果；Ok为对应时刻的观测值。

2 结果与分析

2.1 模型预热轮数与抽样次数M的确定

模型预热以相邻2轮的平均相对误差达0.001为结束条件，判断标准见式(4)。2站点的预热过程类似，此处仅展示禹城站结果（见图1），图1中结果均为各轮预热后模型的最终结果。

图1 禹城站预热过程Fig.1 Spin-up process of Yucheng station

式中：Rk为第k 个模型运行步长时刻的输出结果；K 为模型输出变量的时间维度；n为变量达到稳定所需预热轮数。

不同输出变量达到稳定所需预热时间的结果如表3所示。结果表明，土壤含水率（SM1和SM2）的模拟结果达到稳定所需的时间最短，热通量（LE、H和土壤热通量G）相较需要更长的时间，其中显热通量的预热结果不稳定，预热多轮后其相对误差成周期性循环，取进入循环前一年作为最终所需预热轮数；NEE 模拟稳定所需时间需要100 a 以上，与其值较小（10-4量级），达到数值稳定更为困难不无关系。此结果与密西西比流域的研究结果存在差异，CAI 等[5]认为土壤水预热需要8 a，而显热与潜热通量预热的时间更短，大约需要4 a。

图2 为不同抽样次数的输出变量敏感参数分析结果。结果表明：对于不同抽样次数而言，随着抽样次数的减少，大部分参数的全局敏感度大小发生变化，尤其是植被相关参数（如QE25、SLA、LTOVRC 等）变化较为显著，全局敏感度较高的参数变化较明显。尽管如此，无论抽样次数如何变化，对输出变量不敏感的参数依旧只有非常小的全局敏感度。因此，在不考虑敏感性参数精确排序的条件下，若只为确定哪些参数更重要、更为敏感，拉丁超立方抽样法（LHS）的抽样次数M 可以取参数数量N 的5 倍左右，这与HUO 等[7]的研究结论具有一致性。下文中M 取为250 以节约计算成本。

图2 不同抽样次数的输出变量参数敏感参数分析Fig.2 Parameter sensitivity analysis of output variables for different sampling number

2.2 站点参数敏感性结果分析

为了对参数的敏感性进行区分，简单地将全局敏感度高于0.10 的参数划分为高敏感性参数，将全局敏感度高于0.05且不高于0.10 的参数划分为中敏感性参数，禹城站与栾城站的敏感性参数划分结果见表4。

表4 禹城站与栾城站敏感性参数汇总（btr = 2 : CLM）Tab.4 Summary of sensitivity parameters for two stations（btr = 2 : CLM）

综合2 个站点的结果，由于LE 与H、SM1 与SM2 涉及的物理过程一致，因此其敏感性参数也大致相同。对LE 与H 来说，模拟过程不仅包含水量平衡还涉及到冠层能量平衡过程，因此它们对植被动态生长过程中的植物器官周转率LTOVRC及25 ℃时叶绿素吸收的光量子效率QE25也很敏感。LE、H与土壤有关的敏感参数主要为BEXP、SMCMAX、SMCREF、SMCWLT，这些参数都与植被作用下的土壤水状态密切相关。故而它们也是SM1 的敏感参数。与SM1 不同的是，由于SM2所处位置更深，土壤水与植被的相互作用更小，因此对凋萎系数SMCWLT 及田间持水量SMCREF 这类与植物相关的土壤参数敏感性较弱，而对土壤孔径分布指数BEXP 及土壤孔隙率SMCMAX这类土壤物理性质参数更敏感。

NEE 的敏感性参数为微生物呼吸参数MRP 和叶组织转化速率LTOVRC，但NEE 对MRP 的敏感性显著高于其他敏感参数，这一方面可以证实农田生态系统的固碳潜力主要集中在农田土壤[17]，另一方面与模型对农田生态系统的作物碳同化过程模拟能力较弱有关，Noah-MP 模型的模拟结果低估了冬小麦与夏玉米生长旺季（4 月和8 月）时农田生态系统的碳汇作用（见图3）。同时，LTOVRC 是一个与叶片物质量更新周转进入土壤有机碳有关的参数，这同样能说明在模型模拟农田生态系统NEE 过程中，农田土壤碳库的变化过程占主导，而与光合作用碳同化过程有关的作物碳库变化过程较弱。与光合作用相关的VCMX25、QE25 及SLA，虽然也是重要参数，但NEE 对它们敏感性较低。

图3 禹城站与栾城站NEE模拟值与真实值的比较Fig.3 Comparison of simulated and observed NEE for two stations

整体上，2个站点的敏感参数基本相同（见图4）。结合所有输出变量来看，2 个站点的区别体现在土壤孔隙率SMCMAX、田间持水量SMCREF、凋萎系数SMCWLT 和土壤孔径分布指数BEXP 这4 个参数上。相比禹城站，栾城站模拟的LE、H 对SMCMAX 的变化更不敏感，而SM1、SM2 受SMCREF 和SMCWLT 的影响更大。SMCMAX 用于确定土壤水达到静水平衡时的地下水位，其值为划分土壤水与地下水的阈值，体现了土壤水与地下水间的转换关系。禹城站SMCMAX 更为敏感，说明禹城站的地下水位更浅，地下水对土壤水分状态的影响较大。站点资料显示禹城站地下水位埋深一般为1.5～4.0 m，栾城站地下水位埋深大部分在15 m以下。SMCWLT 用来刻画水分胁迫对植物生长、器官凋亡等生理过程的影响，SMCREF则决定了降水与灌溉水的下渗与土壤的疏干过程，BEXP 为土壤水分特征方程的重要参数，与土壤含水状态密切相关。栾城站的浅层土壤水（SM1、SM2） 对SMCWLT、SMCREF、BEXP 更为敏感，说明栾城站的气候更为干旱（禹城站该轮作年干旱指数为1.929，栾城站为2.540），土壤含水率受降水与灌溉的影响大，作物生长更可能受到水分胁迫，由气象数据根据FAO 推荐的Penman-Monteith 方法计算的日参考作物需水量与降水见图5。

图4 禹城站与栾城站各输出变量参数敏感参数分析（btr = 2 : CLM）Fig.4 Parameter sensitivity analysis of each output variable for two stations

图5 禹城站与栾城站的日参考作物需水量与降水Fig.5 Daily reference crop evapotranspiration and precipitation for two stations

2.3 气孔阻抗的土壤水分限制因子参数化方案的比较

在植被覆盖情况下，气孔阻抗不仅决定了植物蒸腾过程中水汽进入大气的能力，而且也决定了植物的光合作用中二氧化碳进入叶肉组织的能力。植物的气孔开闭受到土壤水分状态的调控，在土壤含水率较低时，植物能感应根区土壤水

分状态并对气孔进行调节，防止自身过分失水而死亡。在Noah-MP 模型中，土壤水分对气孔开闭影响的描述存在于气孔阻抗的土壤水分限制因子参数化方案之中。2.2 节中选用的是基于土壤水势的CLM 方案，本节将结合2.2 节的计算结果，比较CLM 方案与基于土壤含水率的Noah 方案在参数敏感性上的异同。与2.2节的划分原则一致，禹城站与栾城站的中/高敏感参数见表5。考虑参数抽样引起的变异性，在比较不同气孔阻力参数化方案时，与比较站点间差异一样，重点关注全局敏感度的相对大小变化，分析那些相对变化明显的参数。

表5 禹城站与栾城站敏感性参数汇总（btr = 1 : Noah）Tab.5 Summary of sensitivity parameters for two stations（btr = 1 : Noah）

结合2 个站点来看（图6、图7），使用不同气孔阻力参数化方案对SM1、SM2、NPP 的参数敏感性几乎无影响。而将参数化方案从CLM 更换为Noah 2 个站点都出现了以下现象：①对LE 来 说LTOVRC、QE25的敏感性降低，SMCREF、SMCWLT 的敏感度升高；②对H 来说QE25 的敏感性降低。由于控制气孔阻力的土壤水分因子参数化方案的本质区别在于Noah 方案基于土壤含水率而CLM 方案基于水势，在Noah 方案中气孔阻抗由SMCREF 和SMCWLT 直接决定，因此Noah 方案中LE 对这2 个参数更敏感。模型中的水碳耦合的关键枢纽在于叶片的气孔，气孔开放程度不仅影响植物蒸腾还影响光合过程，进而间接影响作物的生长。LTOVRC、QE25 的敏感性降低能在一定程度上反映出LE 对土壤水势的响应灵敏度高于土壤含水率。H 与水通量的联系没有LE 那么紧密，因此与土壤水直接相关的SMCREF、SMCWLT的敏感度变化不明显。而H 对QE25 的敏感性降低但LTOVRC 的敏感度变化不明显，这反映出气孔开闭对H 的影响主要是通过影响植物的光合过程实现。对比图6 与图7 可以看出，在气候更为干旱的栾城站，不同参数化方案下BEXP 的敏感度发生了非常明显的变化，原因在于BEXP直接参与CLM方案的水势计算过程，因此在气候干旱地区使用该方案时需要重视对BEXP的校准。

图6 禹城站不同参数化方案的各输出变量参数敏感参数分析Fig.6 Parameter sensitivity analysis of each output variable for different parameterization schemes at Yucheng station

图7 栾城站不同参数化方案的各输出变量参数敏感参数分析Fig.7 Parameter sensitivity analysis of each output variable for different parameterization schemes at Luancheng station

3 讨 论

华北平原主要的土地利用类型为耕地，是一个大型的农田生态系统，水热碳通量的变化受到种植、管理等人类活动的影响。与其他植被覆盖类型相比，陆面模式在模拟农田时存在着巨大的挑战。为了准确刻画华北平原农田生态系统中的水、能量及CO2通量变化过程，降低参数不确定性对陆面模式的影响，需要开展模型参数敏感性分析。

本研究在ARSENAULT 等[6]的基础上对用于敏感性分析的单站点模型进行修改，完善预热前后浅层土壤中的速效碳库FASTCP、深层土壤中的迟效碳库STBLCP 等模型状态变量的继承关系，代码修改前后栾城站的敏感性分析结果见图8。结果显示，模型代码修改前后，输出变量LE、H、SM1 与SM2对不同参数的敏感程度基本一致，NEE 的敏感性参数发生巨大变化。模型代码修改后，LE、H 对土壤参数的敏感性得以凸显，揭示了热通量与土壤水的潜在联系，参数BEXP 成为最敏感参数，与HUO 等[7]的研究结果一致。BEXP 是土壤的形状参数，与土壤的孔隙结构有关，同时其作为土壤水分特征方程的重要参数，能够反映土壤的持水性，与陆气间水分交换息息相关。对于NEE 来说，模型代码修改后，敏感性参数不再是植物生理相关参数，微生物呼吸参数MRP 的全局敏感性指标显著高于其他参数，这恰能体现农田生态系统碳汇功能主要依赖土壤固碳而不是农作物生物量固碳。同时，本研究发现为使净生态系统碳交换量（NEE）达到稳定值需要经历100 a 以上的预热过程。修改代码并对模型进行充分预热后，为土壤微生物呼吸提供能量的浅层土壤中的速效碳库FASTCP不断积累，这也导致MRP 的敏感性增强。值得注意的是，NEE 预热稳定后FASTCP 的值为102量(t/hm2)，远高于第2 次全国土壤普查数据中农田表层土壤有机碳密度（26.6～32.5 t/hm2）[17]，这是Noah-MP 陆面模式对NEE 等碳循环物理量模拟能力欠缺导致的，开发者应重视对碳循环模块的优化。在对比2种气孔阻抗的土壤水分限制因子参数化方案对模型敏感性的影响时发现，基于土壤水势的CLM 方案能更好地将干旱地区植物根区的缺水状态反映到LE 的变化上；而YANG 等[18]也得出了类似的结论，认为基于土壤含水率的Noah 方案在降水充沛的西双版纳站能获得最接近观测值的模拟结果。ARSENAULT[6]等同样发现在半干旱的草原与农田生态站CLM方案对BEXP 的敏感度高于Noah 方案，并认为PSISAT 在不同方案间的敏感度差异也源自同样的机理。本研究并未发现PSISAT 在不同方案间的敏感度差异，这可能是因为选择的下垫面单一，只聚焦于华北平原生态系统。

图8 代码修改前后栾城站参数敏感性分析结果Fig.8 Parameter sensitivity analysis results at Luancheng station before and after code modification

本研究揭示了Noah-MP 模型模拟华北平原农田生态系统水热碳通量及土壤含水率时的敏感参数，并对比了不同气孔阻抗的土壤水分限制因子参数化方案下参数敏感度，对参数优化与参数化方案选择有着指导意义。但本研究仍存在一些不足，例如，只考虑用户自定义参数的敏感性，而未考虑模型嵌入参数，CUNTZ 等[19]发现对潜热通量、总径流及各径流分量而言，许多模型嵌入参数甚至比用户自定义参数更敏感，因此未来可以将参数敏感性分析的范围拓展至模型嵌入参数。

4 结 论

（1）若不需要对参数按敏感性强弱进行严格排序，在使用拉丁超立方抽样法（LHS）筛选模型的敏感性参数时，为节约计算成本，抽样次数可以取为参数的5倍左右。

（2）用Noah-MP 模型模拟华北平原农田生态系统典型站点时，LE、H 的敏感参数有：LTOVRC、QE25、BEXP、SMCMAX、SMCREF、SMCWLT。其中，BEXP、SMCMAX、SMCREF、SMCWLT 也是SM1 和SM2 的敏感参数。NEE 对MRP、LTOVRC 的敏感性显著高于VCMX25、QE25 及SLA 等与作物光合作用相关的参数，体现农田生态系统碳汇功能主要依赖农田土壤固碳。

（3）在气候干旱的地区需要重视对土壤参数SMCMAX、SMCREF、SMCWLT、BEXP 的校准，尤其在使用CLM 气孔阻抗的土壤水分限制因子参数化方案时，BEXP 的校准需要特别关注。

致谢：感谢禹城农业综合试验站与栾城农业生态系统试验站为本研究提供的数据支持，感谢Arsenault Kristi R.在GitHub 平台贡献的开源代码，使本研究能在此基础上提出改进。