基于VIP 模型评估不同全局敏感性分析方法有效性及效率

何力鸿，莫兴国，刘苏峡

（1. 中国科学院地理科学与资源研究所陆地水循环及地表过程重点实验室，北京 100101；2. 中国科学院大学，北京 100049；3. 中国科学院大学资环学院/中丹学院，北京 100049）

0 引 言

基于动力学过程的生态系统模型具有结构复杂、输入参数多和参数空间变异性强等特点，使得参数校验成为模型研究和应用的难点[1]。而参数敏感性分析可有效识别模型敏感参数，通过固定不敏感参数，可有效降低运算成本[2]。

参数敏感性分析可分为局部和全局敏感性分析。其中全局敏感性分析方法充分考虑不同参数及其相互作用对模型模拟结果的影响，在非线性过程模型的参数校准中得到广泛应用[3-5]。全局敏感性分析方法可分为定性和定量分析，其中定性敏感性方法研究不同参数对目标函数影响的相对大小，如 Morris 法[6]、多元自适应回归样条法（Multivariate Adaptive Regression Splines，MARS)[7]、Delta Test 法[8]和 Sum of Trees 法[9]等；而定量敏感性方法可获得不同参数对目标函数的贡献率，如 Sobol'法[10]、McKay 法[11]和扩展傅里叶幅度检验法 （Extended Fourier Amplitude Sensitivity Test，EFAST)[12]等。近年来，国内外学者对参数敏感性做了诸多分析，如梁浩等[13]通过Morris和Sobol'法分析WHCNS 模型中与土壤硝态氮和含水率相关的敏感参数；邢会敏等[14]采 EFAST 法对AquaCrop 模型中的作物参数进行分析。DeJonge 等[15]使用 Morris 和 Sobol'方法分析了不同灌溉处理下CERES-Maize 模型中作物参数对叶面积指数、蒸散和产量的影响，并基于敏感度分析结果提出了参数校准方案。

在参数分析过程中，合适的试验设计可减少抽样次数，降低运算成本。目前试验设计方法主要有蒙特卡洛（Monte Carlo，MC）法[16]、正交试验设计（Orthogonal Experimental Design，OED）[17]和拉丁超立方（Latin Hypercube，LH）抽样[18]等。其中，MC 通过参数空间上的概率密度函数随机生成样本点，其效率取决于随机数的属性。而随机数抽样容易产生聚类，出现空白抽样空间，随机添加新抽样点不一定能填补空白抽样空间[2]。有限数量的抽样点的正交质量取决于点分布的均匀性，为了保证所有的抽样区域都能被抽样点覆盖，McKay 等[11]提出了LH 抽样，该方法是分层抽样的一种形式，可减少通过蒙特卡洛法估算时产生的误差。后来研究对LH 抽样进行了改进，如Tang[19]通过正交阵列（Orthogonal Array，OA）构建了OALH（Orthogonal Array-Based Latin Hypercubes），与LH 相比，OALH 更适用于计算机试验和数值积分。除上述抽样方法外，还有针对特定敏感性方法设计的抽样方式，如 Morris 抽样[6]和Sobol'抽样[10]等。

以往研究多基于一种或几种敏感性分析方法筛选敏感性参数，对比分析不同敏感性方法有效性及效率的研究不够深入。本研究基于VIP（Vegetation Interface Processes）模型探讨不同敏感性方法和抽样方式在参数敏感性分析过程中的有效性及其效率，VIP 模型综合考虑了作物生长、碳氮循环过程、土壤水分运动等过程，经过多年发展，模型已用于华北平原单点及区域作物产量、蒸散量及水分利用效率的模拟[20]。本研究以华北平原土壤硝态氮模拟为例，对比分析了PSUADE（Problem Solving environment for Uncertainty Analysis and Design Exploration）提供的8 种敏感性分析方法（包括6 种定性和2 种定量敏感性分析方法）在筛选农田土壤硝态氮相关敏感参数时的有效性及其效率。

1 材料与方法

1.1 VIP 模型及参数

VIP 模型综合考虑了植物生理过程、碳氮循环过程、土壤水分运动等过程，能够合理模拟单点和区域作物生长和水量平衡过程。植被动态变化采用干物质累积与分配方程来描述、土壤热量传输采用热扩散方程、土壤水量平衡方程遵循达西定律[20]。其中，VIP 模型土壤氮循环模块考虑了土壤有机氮库之间的周转、矿化过程、无机氮库之间的转化过程。本研究选取了与土壤有机质周转、矿化、硝化过程、反硝化过程和氨挥发过程相关的参数（表1）来评估不同敏感性分析方法。

表1 VIP 模型中氮循环过程相关参数Table 1 Parameters related to nitrogen cycle process in VIP(Vegetation Interface Processes) model

1.2 试验站点及数据

研究样地位于华北平原 (31°14'～40°25'N，112°48'～122°45'E)，面积约33 万km2，其中农业用地占70%。属半干旱半湿润暖温带大陆性季风气候，夏季降雨量最多[21]。中部地区的年平均气温在10～15 ℃之间，日照时间约为2 800 h。土壤为黄河沉积物发育而成的轻壤质潮土，耕层容重为1.48 g/m3，有机质质量分数为11 g/kg。该地区玉米于每年6 月初播种，9 月中下旬收获。本研究选择了3 个站点（封丘，禹城和栾城农业生态站）的数据对VIP模型进行参数敏感性分析。

1.2.1 驱动数据

气象数据、土壤及作物参数等数据用于驱动 VIP 模型，其中气象数据包括最大温度、最小温度、平均温度、降水量、日照时数、相对湿度、风速。封丘，禹城和栾城站日尺度气象数据从国家生态科学数据中心[21]获得。另一类驱动数据是作物参数与土壤参数，其中涉及的作物参数主要有：有效积温参数、作物不同时期含氮量；土壤参数主要有不同有机质分解速率、碳氮比、氨挥发速率、反硝化速率的相关参数，由文献资料获得[22-23]。封丘、禹城和栾城农业生态站田间管理措施数据从文献资料[24-25]和国家生态科学数据中心[21]获得。表 2 中列出了播种日期，收获日期，氮肥和灌溉的详细信息。考虑到研究目的、数据的质量等因素，本研究封丘站、禹城站和栾城站模拟时间段分别为2008年6月—2009年9月，2001 年 6 月—2002 年 9 月和 2004 年 6 月—2005 年 9 月，时间步长定为1 d。

表2 3 个农业生态站玉米生长期及田间管理数据Table 2 Growing season and field management data of maize at three agro-ecological stations

1.2.2 验证数据

本研究选定 3 个农业生态站（封丘、禹城和栾城）样地实测土壤硝态氮数据用于验证 VIP 模型模拟结果。在玉米生长季，以20 cm 的距离分层收集土壤样品，并在分层中的 2 个点随机采集每个样区。在每季作物收获后采集土样，用50 mL 2 mol/L KCL 提取10 g 每个土壤样品，采用双波长（220 nm，275 nm）法测定[26]。

1.3 敏感性分析方法

在以往农业生态系统模型参数敏感性分析研究中，多数研究基于一种或几种敏感性分析方法筛选敏感性参数，而对比分析不同敏感性方法有效性及效率研究较少。本研究选取了在生态模型中广泛应用的 Spearman秩相关系数（Spearman’s correlation coefficient，SPEA）法、Morris 法、多元自适应回归样条（Multivariate Adaptive Regression Splines，MARS)、Delta Test（DT)法、Sum of Trees（SOT）法、Gaussian Process（GP）法、McKay 法和Sobol'法，对比分析这8 种全局敏感性分析方法在参数敏感性分析中的有效性及效率。8 种全局敏感性分析方法由 PSUADE（Problem Solving environment for Uncertainty Analysis and Design Exploration）软件提供，PSUADE 为复杂系统模型提供不确定性分析的集成设计和分析环境，可通过 https://computing.llnl.gov/projects/psuade-uncertainty-quantification获得。在参数敏感性分析过程中，首先指定抽样方式，样本量大小以及参数取值范围，然后生成样本组合。本研究使用抽样方式有蒙特卡罗（Monte Carlo，MC）、METIS、Morris、拉丁超立方（Latin Hypercube，LH）、LPTAU、正交阵列 （Orthogonal Array，OA)和基于拉丁超立方的正交阵列（Orthogonal Array-Based Latin Hypercubes，OALH）抽样方式。通过脚本将 VIP 模型可执行文件（vip.exe）与PSUADE 进行链接，并将运行结果读写到输出文件，利用PSUADE 提供的统计方法分析不同参数组合与目标函数之间关系。

在评估敏感性分析方法有效性时，将敏感性分析结果进行归一化，通过对比不同敏感性分析结果及关键参数的筛选来判断敏感性方法的有效性，如潜在硝化速率（knit）对土壤硝态氮含量影响较大，若筛选结果判定其为敏感参数，则证明该敏感性分析方法有效。在评估敏感性方法效率时，通过筛选敏感参数所需样本量大小来衡量不同敏感性方法效率，样本量越大，效率越低，反之效率越高。在评估Morris 敏感性方法效率时，样本量是n+1（n是参数个数）的整数倍，为探究筛选敏感参数所需最佳样本量，设置样本量分别为85, 170,…, 3400，每次增加 85 直至最大样本量；Sobol'敏感性分析方法筛选敏感性参数过程中，样本量是n+2（n是参数个数）的整数倍，本研究设置样本量分别为90, 180,…, 3420，每次增加 90 直至最大样本量；在评估 DT、MARS、Mckay和SOT 敏感性方法效率时，通过设置不同抽样方式（MC，LH，OA，OALH，LPTAU 和METIS）和不同样本量梯度分析有效敏感性分析方法的效率，其中 MC，LH，LPTAU 和 METIS 抽样方式的样本量设置为 85, 170,…,3400，样本量每次增加85 直至最大样本量；OA 和OALH抽样法的样本量为m·pt，本研究给定m=1 和t= 2，样本量由素数p（17, 19, 23,…, 59）确定，样本量分别为289,361, 529,…, 3481。样本量从小依次增加到最大值（最大样本量分析结果作为基准），对比不同样本量下敏感性分析结果与最大样本量的结果，确定筛选出敏感参数所需的最少样本量。

2 结果与分析

2.1 敏感性分析方法有效性

2.1.1 定性敏感性分析方法有效性

将 3 个农业生态站定性敏感性方法（SPEA、DT、GP、MARS、Morris 和SOT）的分析结果进行加权平均，为了方便统一对比，对结果再进行归一化处理，本研究将敏感性分析归一化结果大于0.2 的参数作为敏感参数。对VIP 模型中16 个氮循环参数进行排序（图1），结果表明，DT、MARS、Morris 和SOT 法筛选出与土壤硝态氮含量相关的敏感参数有：knit、kur、mnit、rmb、rsh、kmb、ksl和kml。其中，knit和kur在上述4 种敏感性分析方法筛选结果中表现出高敏感性，而SPEA 和GP 法筛选结果表明knit和kur为不敏感性参数，与其他定性敏感性分析结果不一致。此外，SPEA 法也无法有效识别kmb，GP 法无法识别ksl，而参数kmb和ksl在 DT、MARS、Morris 和 SOT法的结果中都表现为敏感参数。

2.1.2 定量敏感性分析方法有效性

在 McKay 和 Sobol'定量敏感性分析方法的有效性（图 2）验证过程中，2 种方法样本量都设置为3 400。为了方便统一对比，将 3 个农业生态站定量敏感性分析结果进行归一化处理。结果表明，McKay 法和 Sobol'法都能有效筛选出 3 个农业生态站土壤硝态氮相关的敏感参数，敏感参数分别为knit、kur、mnit、ksl、kml、ksh、kmb、rmb和rsh，且筛选出的最敏感参数都是knit，这与定性敏感性分析结果一致。

图1 不同定性敏感性分析方法的参数敏感性Fig.1 Parameter sensitivity of different qualitative sensitivity analysis methods

2.2 敏感性分析方法效率

由于SPEA 和GP 法敏感性分析结果与其他敏感性分析方法结果不一致，未能将最敏感参数knit有效筛选出来，未进一步分析SPEA 和GP 法效率。本研究通过6 种不同的抽样方法（MC，LH，OA，OALH，LPTAU 和METIS）来明确适用于MARS、DT、SOT 和Mckay 敏感性分析法的抽样方式及最小样本量。在对比不同敏感性分析方法效率时，将每种抽样方式最大样本量筛选结果作为敏感性方法有效性基准，当连续 3 个样本量的敏感参数与最大样本量筛选结果相同或接近时，则认为连续3 个样本量中最小样本量足以用于抽样分析。对不同抽样方式及样本量的分析结果进行统计（表3），结果发现，MARS敏感性分析过程中，通过 OA 和 OALH 抽样时，所需样本量仅为361，这表明OA 和OALH 抽样方式最适合MARS 法；在利用DT 和SOT 敏感性分析中发现，MC最适合DT 和SOT 敏感性分析方法，所需最小样本量分别为425 和510；与MARS 和DT 敏感性分析方法相比，SOT 在筛选敏感参数过程中所需样本量更多。在McKay 筛选敏感参数的过程中发现，与其他抽样方式相比，OALH 抽样需要的样本点最小（样本量为529），这表明，OALH 抽样方式最适应于McKay 敏感性分析方法。

图2 3 个农业生态站McKay 和Sobol'敏感性分析方法的参数敏感性Fig.2 Parameter sensitivity of Mckay and Sobol' methods at three agro-ecological stations .

表3 3 个农业生态站不同抽样方式下MARS、DT、SOT 和Mckay 敏感性分析方法所需最小样本量Table 3 Minimum sample size required for MARS, DT, SOT and Mckay sensitivity analysis methods under different sampling methods at three agro-ecological stations

在分析Morris 法和Sobol'法筛选敏感参数效率时，设置不同样本量来确定筛选华北平原 3 个农业生态站（封丘、禹城和栾城）土壤硝态氮敏感参数所需的最小样本量。不同参数对应均值μ在4 kg/ hm2处会被明显区分开，故本研究认为均值μ≥ 4 kg/ hm2的参数为敏感参数。结果表明，样本量为 85 时，封丘、禹城和栾城站筛选出的敏感参数为knit和kur；随着样本量增加，mnit、ksl、kml、ksh、kmb、rmb和rsh被筛选出来；样本量大于 340筛选出的敏感参数基本一致，这表明样本量取340 时即能有效筛选 VIP 模型中土壤氮循环敏感参数。在分析Sobol'法筛选敏感参数效率时，在分析Sobol'法筛选敏感参数效率时，当参数 Sobol'一阶指数大于 0.4 时，认为该参数为敏感参数。通过对比不同样本量下敏感性分析结果发现（表 4），除了禹城生态站，当样本量为 720时，封丘和栾城生态站的kmb、ksh都不能被有效筛选出来；而样本量大于810 的敏感性分析结果基本一致，筛选出的敏感参数分别为knit、kur、mnit、ksl、kml、ksh、kmb、rmb和rsh，这表明该研究中样本量选取 810 就能有效筛选敏感参数。

表4 3 个农业生态站Morris 和Sobol'敏感性分析方法所需最小样本量Table 4 Minimum sample size required for Morris and Sobol'sensitivity analysis methods at three agro-ecological stations

3 讨 论

3.1 不同敏感性分析方法有效性

8 种全局敏感性方法运用于土壤氮循环参数进行了全局敏感性分析。从全局敏感性分析结果来看，SPEA 和GP 法敏感参数筛选结果与其他方法的结果不同，主要由于 SPEA 指数是一种秩变换的统计量，适用于参数和模型响应之间的线性关系，而大多数生态系统模型中参数与模拟结果呈非线性关系，如Duan 等[27]分析了SIXPAR模型模拟结果对参数的响应，尽管 SIXPAR 模型结构简单，但结果表明不同参数之间表现出强烈的非线性交互作用。Gan 等[2]通过SAC-SMA 模型对比不同敏感性分析方法时发现，SPEA 和GP 法无法有效筛选出敏感参数，这与本研究结果类似。此外，与SOT 法相比，通过MARS和DT 法筛选出的敏感参数结果接近，且部分边缘敏感参数敏感性高于SOT 法得到的结果，由于边缘敏感参数与不敏感参数之间存在较大差异，可以推断 MARS 和 DT法可能更有效识别不敏感参数；与MARS 和DT 法相比，SOT 法筛查所采用的不同抽样方式所需的样本量差异更大。这与Gan 等[2]研究结果类似。无论采用何种抽样方式，SOT 法筛选比MARS 和DT 法筛选需要更多的样本量来识别不敏感的参数，这支持了MARS 和DT 法擅长识别不敏感的参数，而SOT 法擅长识别高度敏感的参数的推断。此外，与上述定性敏感性分析法相比，McKay 和Sobol′定量敏感性分析法都能有效识别敏感参数，但这2 种方法在敏感性分析过程中运算时间较长，所需计算成本较高。

3.2 不同敏感性分析方法效率

本研究为了确定有效筛选的最小样本量，在使用Morris 方法时，当样本量低于 340 时，部分敏感参数未能被筛选出来，这表明样本量低于 340 时不能有效筛选VIP 模型中氮循环敏感参数，这与以往研究结果类似[28]。通过对比不同样本量下Sobol'一阶指数结果发现，当样本量小于810 时，rmb和rsh不能被有效筛选出来。与Sobol'法相比，Morris 方法所需样本量较小，对于复杂碳氮模型，随着样本量的增加，所需运算成本也增加，而Morris法在样本量较少的情况下就能筛选出较敏感参数，所以在初步筛选复杂模型的敏感参数时推荐使用 Morris 方法。McKay 和 Sobol′方法的敏感性指标丰富，可以从多角度考虑，得到可靠的结果，但McKay 和Sobol′方法需要基于较大样本量才能得到可靠的定量敏感性分析结果，不适用于参数较多的复杂动力系统模式的参数筛选。在实际的应用中，要根据研究目的选择最恰当的方法。如需定量分析参数对模拟结果影响，并不考虑模型运算成本问题，推荐使用McKay 和Sobol′方法。

4 结 论

本研究以VIP（Vegetation Interface Processes）模型模拟华北平原 3 个农业生态站（封丘、禹城和栾城）土壤硝态氮为例，评估不同敏感性分析方法有效性及其效率。通过不同定性敏感性分析结果发现，SPEA（Spearman’s correlation coefficient）和 GP（Gaussian Process）法筛选结果与 DT（Delta Test)、MARS（Multivariate Adaptive Regression Splines）、Morris 和 SOT（Sum of Trees）法的不一致；与MARS、DT 和SOT 定性敏感性分析方法相比，Morris 在筛选敏感参数所需样本量最少；蒙特卡罗（Monte Carlo，MC），正交阵列（Orthogonal Array，OA）和基于拉丁超立方体的正交阵列（Orthogonal Array-Based Latin Hypercubes，OALH）抽样方式都适用于MARS、DT 和SOT 敏感性分析方法。OA 和OALH 抽样方式最适合MARS 法；MC 最适合DT和SOT 敏感性分析方法；OALH 最适应于McKay 敏感性分析方法。与MARS 和DT 敏感性分析方法相比，SOT在筛选敏感参数过程中所需样本量更多。

在筛选敏感性参数过程中，定性敏感性分析方法比定量敏感性方法所需样本量少，但无法对敏感参数进行定量描述。对参数较多的复杂系统模型，可首先使用定性敏感性分析方法进行初步敏感性参数筛选，以较低的计算成本剔除不敏感的参数，然后采用定量敏感性分析方法对系统模型中的敏感参数进行定量分析。