全 文 :生 态 模 型 的 灵 敏 度 分 析 3
徐崇刚1 ,2 3 3 胡远满1 常 禹1 姜 艳3 李秀珍1 布仁仓1 贺红士1
(1 中国科学院沈阳应用生态研究所 ,沈阳 110016 ;2 中国科学院研究生院 ,北京 100039 ;
3 华盛顿大学森林资源学院 ,西雅图 9819522100)
【摘要】 灵敏度分析用于定性或定量地评价模型参数误差对模型结果产生的影响 ,是模型参数化过程和
模型校正过程中的有用工具 ,具有重要的生态学意义. 灵敏度分析包括局部灵敏度分析和全局灵敏度分
析.局部灵敏度分析只检验单个参数的变化对模型结果的影响程度 ;全局灵敏度分析则检验多个参数的变
化对模型运行结果总的影响 ,并分析每一个参数及其参数之间相互作用对模型结果的影响. 目前 ,在对生
态模型的灵敏度分析中 ,越来越倾向于使用全局灵敏度分析的方法. 但国内仍多采用局部灵敏度分析方
法 ,很少采用全局灵敏度分析方法. 文中详细论述了局部灵敏分析和全局灵敏度分析的主要方法 (一次变
换法、多元回归法、Morris 法、Sobol’法、傅里叶幅度灵敏度检验法和傅里叶幅度灵敏度检验扩展法) ,希望
能为国内生态模型的发展提供一个比较完善的灵敏度分析方法库. 结合国内外的灵敏度分析发展现状 ,指
出联合灵敏度研究、灵敏度共性研究及空间直观景观模型的灵敏度分析将为生态模型灵敏度分析研究中
的热点和难点.
关键词 生态模型 灵敏度分析 局部灵敏度分析 全局灵敏度分析 空间直观景观模型
文章编号 1001 - 9332 (2004) 06 - 1056 - 07 中图分类号 Q14 文献标识码 A
Sensitivity analysis in ecological modeling. XU Chonggang1 ,2 , HU Yuanman1 , CHAN G Yu1 , J IAN G Yan3 ,
L I Xiuzhen1 ,BU Renchang1 , HE Hongshi1 ( 1 Institute of A pplied Ecology , Chinese Academy of Sciences ,
S henyang 110016 , China ;2 Graduate School of Chinese Academy of Sciences , Beijing 100039 , China ;3 College
of Forest Resources , U niversity of W ashington , Seatle 9819522100) . 2Chin. J . A ppl . Ecol . ,2004 ,15 (6) :1056
~1062.
Sensitivity analysis is used to qualitatively or quantitatively apportion the variation of model output to different
source of variation. It is a very useful tool in model parameterization and calibration ,and has important ecological
significance by identifying the governing factors for a certain ecological process simulated. There are two schools
of sensitivity analysis , local sensitivity analysis and global sensitivity analysis. The former examines the local re2
sponse of the output ( s) by varying input parameters one at a time ,holding other parameters to a central value ;
and the latter examines the global response (averaged over the variation of all the parameters) of model output
( s) by exploring a finite (or even an infinite) region. Since it is very easy to conduct local sensitivity analysis ,it
is very popular in ecological models. However ,local sensitivity analysis is not computationally effective ,because it
can only get the sensitivity of a single parameter at a time. It can not take into consideration the effect of interac2
tion of different parameters. Additionally ,the value of other parameters will affect the sensitivity of the parameter
specified. In view of this ,global sensitivity analysis is increasingly preferred to local sensitivity in recent years.
However ,for most of the ecological modeling study published in Chinese ,only local sensitivity analysis is conduct2
ed. To provide a toolbox of alternative sensitivity analysis algorithms for ecological model development in future
study ,we reviewed the main methods of both local sensitivity analysis and global sensitivity analysis ,including
one at a time method ,multivariate regression ,Morris’method , Sobol’s method , Fourier Amplitude Sensitivity
Analysis ,and Extended Fourier Amplitude Sensitivity Analysis. Based on the state2of2the2art research on sensitiv2
ity analysis ,the sensitivity of the interaction between two or more than two model parameters ,the sensitivity of
common model parameters in a set of models ,and the sensitivity analysis in spatially explicit landscape model sim2
ulation were identified as the key areas and difficulties of future study on sensitivity analysis in ecological modeling.
Key words Ecological model , Sensitivity analysis , Local sensitivity analysis , Global sensitivity analysis , Spa2
tially explicit landscape model. 3 中国科学院引进国外杰出人才项目 (BR010403) ,国家自然科学基
金项目 ( 40331008 ) 和中国科学院知识创新工程资助项目
(SCXZY0102) .3 3 通讯联系人.
2003 - 03 - 24 收稿 ,2004 - 01 - 04 接受.
1 引 言
生态模型是指对生态现象和生态过程进行模拟的计算
机程序或数学方程 . 在生态模型的发展初期 ,模型的建立者
一般不考虑模型参数的估计误差. 实际上 ,测量误差的存在
和对机理的不完全理解使模型参数具有很大的不确定性. 模
型参数的不确定性必然会使模型运行的结果存在不确定性.
应 用 生 态 学 报 2004 年 6 月 第 15 卷 第 6 期
CHIN ESE JOURNAL OF APPL IED ECOLO GY ,J un. 2004 ,15 (6)∶1056~1062
如果模型结果的不确定性很大 ,模型结果就不能作为可靠的
决策依据. 若要提高模型预测的精度 ,就需要提高模型各参
数的精度 (降低模型各参数的不确定性) . 然而很多生态模型
存在几十到几百个参数 ,要提高每一个参数的精度则很难做
到. 此外 ,由于自然界是一个非常复杂的系统 ,每一个生态过
程都受各种各样的不确定性因素影响 ,所以某些参数的不确
定性是无法降低的. 比如 ,一个种群中每一个个体的出生或
死亡都会受当时的食物、天气状况、捕食者的多少、疾病的传
播及其它各种不确定性因素的影响 ,种群增长率或死亡率的
估计总会存在较大的不确定性. 因此 ,需要通过灵敏度分析
(sensitivity analysis)来评价各个参数的不确定性对模型运行
结果的影响 ,集中人力、物力提高那些对模型结果影响程度
大的参数的精度 ,对于那些对模型结果影响不大的参数 ,只
需选取其经验值.
灵敏度分析也是模型参数校正过程中的一个非常有用
的工具[4 ,16 ,17 ,20 ,29 ,39 ] ,其目的在于确定模型中哪些方面最容
易在系统描述中引进不确定性. 通过灵敏度分析可以确定模
型各参数对输出结果影响的大小 ,在模型校正过程中重点考
虑那些对输出结果影响大的参数 ,对于那些对模型结果几乎
没有影响的参数可以不予考虑 ,这会在很大程度上减小模型
校正的工作量.
灵敏度分析具有很重要的生态学意义. 通过灵敏度分
析 ,可知模型对哪些参数的变化敏感 ,从而可以确定各影响
因子对模型所模拟的生态过程的影响程度. 祝增荣等 [52 ]对
一个白背飞虱种群动态模型的灵敏度分析表明 ,白背飞虱种
群对迁入的时间、量和集中程度变化最敏感 ,说明种群迁入
的时间、量和集中程度是白背飞虱大发生的主导因素. Pacala
等[36 ]对一个空间直观的森林景观模型 ( SORTIE) 进行了灵
敏度分析 ,表明模型模拟结果对初始物种的丰富度在前 300
年比较敏感 ,在 300 年之后不敏感 ,说明物种的初始丰富度
对森林演替的影响能延续 300 年左右.
尽管灵敏度分析的重要性被广泛认同 ,然而文献中的很
多模型并没有进行这项分析 ,主要是因为计算过程比较复
杂 ,而且计算结果很难解释 [28 ] . Tomovic′[49 ]第一次用系统综
合的形式描述了灵敏度分析的数学基础. Waide 和 Web2
ster[50 ]详细讨论了将参数变化作为灵敏度分析的方法 ,用以
评价生态系统模型参数的灵敏度. Recknagel[40 ]用 83 个系数
分析了一个生态系统模型 ,表现出灵敏度分析的复杂性和数
学活力. 此后 , 灵敏度分析被广泛用于各种生态模型
中[1 ,15 ,18 ,28 ,30 ,35 ,37 ,38 ,45 ] .
灵敏度分析定性或定量地评价模型参数不确定性对模
型结果的影响[7 ] . 灵敏度分析包括局部灵敏度分析 (local
sensitivity analysis) 和全局灵敏度分析 (global sensitivity anal2
ysis) [42 ] . 局部灵敏度分析只检验单个参数的变化对模型结
果的影响程度 ,其它参数只取其中心值. 全局灵敏度分析检
验多个参数的变化对模型结果产生的总的影响 ,并分析每一
个参数及参数之间的相互作用对模型结果的影响. 全局灵敏
度分析与局部灵敏度分析的区别在于 :1)每一个参数在一个
有限大甚至无限的范围内变化 ;2)由某个参数变化引起的模
型结果的变化是全局的 ,即模型结果的变化是在所有参数变
化的共同作用下产生的. 全局灵敏度分析又可分为定性全局
灵敏度分析和定量全局灵敏度分析. 定性全局灵敏度分析只
是定性地分析模型各参数的不确定性对模型结果影响的相
对大小 ,也称因子筛选灵敏度分析 [42 ] ,其目的是以较低的计
算代价获取模型各个输入参数灵敏度大小的排序. 定量全局
灵敏度分析则定量地给出各参数的不确定性对模型输出结
果不确定性的贡献率.
局部灵敏度分析的优点在于其可操作性. 但是局部灵敏
度分析也存在以下缺点 :1) 只考虑一个参数的变化量 ,一次
只能分析一个参数 ,在计算上不是很有效 ;2) 由于一次只对
一个参数进行分析 ,所以不能考虑模型参数之间的相互作用
对模型输出结果的影响 ;3) 对某一参数进行分析时 ,其它参
数的不同取值会影响其灵敏度. 因此 ,目前对生态模型的灵
敏度分析开始越来越倾向于使用全局灵敏度分析的方法. 但
是在国内对生态模型的灵敏度分析中 ,一般都采用局部灵敏
度分析的方法[51 ,52 ] ,而未见采用目前迅速发展的全局灵敏
度分析. 本文详细论述了灵敏度分析的主要方法 ,希望能为
国内生态模型的发展提供一个比较完善的灵敏度分析方法
库.并结合灵敏度分析研究发展现状 ,指出目前生态模型灵
敏度分析研究中的热点和难点.
2 局部灵敏度分析
局部灵敏度分析也称一次变化法 ,其特点是只针对一个
参数 ,对其它参数取其中心值 ,评价模型结果在该参数每次
发生变化时的变化量. 有两种变换法 :第一种是因子变化法 ,
如将预分析的参数增加 10 %或减少 10 % ;另一种方法是偏
差变化法 ,如将预分析的参数增加一个标准偏差或减少一个
标准偏差. 通常会采用灵敏度系数作为衡量参数灵敏度的标
准. 最简单的灵敏度系数的形式是 :
S i =
dv
dpi
(1)
其中 , S i 是第 i 个参数的灵敏度 , v 是指所预测的模型的结果
参数 , pi 是第 i 个参数. 当然 ,在实际的应用中会对公式 (1)
作一定的扩充[28 ,51 ] .
人们通常认为 ,如果模型对参数误差不敏感 ,则这个模
型就可应用到不同的环境条件中 ,其实这种结论是不正确
的. 因为真正的参数错误可能不止是 10 %或 20 %的不确定
性 ,也可能是 2 倍或 10 倍的不确定性 ,如果模型参数初始化
不正确的话 ,那么模型对于 10 %的变化是不敏感的.
3 定性全局灵敏度分析
311 多元回归法
多元回归法基于拉丁超几何体取样 (Latin hypercube
sampling) . 拉丁超几何体取样法最先由 Conover 在 1975 年
提出 ,并由 Mc Kay 等[31 ]在 1979 年正式发表. 由于把参数的
累积概率分布分成多个等间隔的区间 ,拉丁超几何体取样法
75016 期 徐崇刚等 :生态模型的灵敏度分析
比普通的随机采样法更有效 [23 ] . 它将每一个参数的概率分
布函数的纵轴 (0~1)分成几个互不重叠的等间隔区间 ,每一
个区间分别对应横轴参数的一个等概率取值区间. 因此 ,概
率分布函数的定义域被分成几个互不重叠的等概率取值区
间.在纵轴上的每一个区间内随机采一个样 ,就对应横轴上
的一个参数值. 如果一个模型具有 n 个参数 ,而每一个参数
的概率分布都分成 m 个区间 ,那么将有 nm 个取样组合. 而
实际上 ,只是取 m 个样. 其做法是 :把 n 个参数的取值排列
成一个 n ×m 的矩阵 ,把矩阵每一列元素的次序随机打乱 ,
就得到了所有参数的 m 个输入. 矩阵的每一行作为各个参
数的输入值 ,运行模型获得模型的模拟结果. 用模拟结果与
各输入参数的多元线性回归系数或偏相关系数表示各参数
的灵敏度大小 ,也可用逐步回归的方法确定各参数的灵敏度
大小. 基于拉丁超几何体取样的多元回归法的详细介绍见
Mc Kay 等[31 ] ,Downing 等[13 ]和 Helton[19 ] . 它已在生态模型
的灵敏度分析中获得了广泛应用. Mcarthy 等[30 ]用一个随机
的状态结构模型 ( stochastic stage2structured model) 对蜜雀
( L ichenostom us melanops cassidix)种群进行种群变异性分析
(population variability analysis) 时 ,多元回归灵敏度分析表
明 ,在 4 个生育力参数中未受干扰配偶的平均生育力对种群
所面临灭绝风险的影响最大 ;Nathan 等[34 ]对一个由风引起
的种子传播机理模型进行多元回归灵敏度分析时 ,用逐步回
归法确定模型各参数的灵敏度大小 ,结果表明水平和垂直的
风速对种子传播的影响远远大于物种的生物学特性 ; Shirley
等[46 ]对一个蛞蝓 ( Deroceras reticulatum) 种群模型进行灵敏
度分析时 ,先根据已有的观察数据确定每一个模型参数的分
布 ,采用拉丁超几何体取样法对每个参数随机取样 ,最后用
多元回归的偏相关系数量化每一个参数的灵敏度. 分析表
明 ,对模型输出结果有显著影响的环境参数为土壤湿度、空
气温度和叶面积指数 ;有显著影响的物种生活史参数为 :个
体生长速度、由环境因素导致的死亡率 (environmental mor2
tality)和成熟个体的死亡率.
312 Morris 法
Morris法是由 Morris 在 1991 提出的[33 ] . 它先将每个参
数的取值范围映射到[0 ,1 ] ,并将其离散化 ,使每个参数只从
{ 0 ,1/ ( p - 1) ,2/ ( p - 1) , . . . ,1} 中取值. 其中 , p 为参数的
取样点的个数. 每一个参数都在 p 个取样点上随机取值 ,获
得向量 X = [ x1 , x2 , ⋯, x i , ⋯⋯xk ] ,其中 k 为参数的个数.
考虑 m ×k ( m = k + 1) 的矩阵 B
0 0 0 ⋯ 0
1 0 0 ⋯ 0
1 1 0 ⋯ 0
1 1 1 ⋯ 0
⋯ ⋯ ⋯ ⋯ ⋯
1 1 1 1 1
及变化量 △ = s/ ( p - 1) ( s 为变化因子) ,显然 ,矩阵 △B 中
相邻两行只有一个参数的取值不同 , 且其变化量为 △. 因
此 ,可以把这两行作为模型的参数输入 ,分别获得模型的输
出结果 y1 和 y2 . 公式5 i ( X) = ( y1 - y2) / Δ (2)
可用于计算模型参数 i 的灵敏度 (为了便于比较 ,在实际计
算中可取其绝对值) . 取所有 k 组相邻行元素作为模型的输
入参数可以获得所有 k 个参数的灵敏度. 由此可见 ,只需要
一次随机取样就可以得到 k 个参数的灵敏度 ,这在计算上是
很有效的. 但矩阵 △B 中各元素的取值并不是随机的 ,在实
际的操作中 ,对矩阵 B 采用一随机过程以保证取值的随机
性. 该随机过程如下 :
1) 设 D 3 为 k 维对角阵 ,对角线元素为 + 1 或 - 1 的概
率相等. 设 J m , k 为 m × k 的单位矩阵 , 矩阵 (1/ 2) [ (2B -
J m , k) D 3 +J m , k ]为 m ×k的矩阵 ,矩阵中每一列元素与 B 中
相对应的元素相等 ,或把 0 替换为 1 ,或把 1 替换为 0 ;
2) 设 X3 为 X的“基值”向量 ,X 中的每一个参数都从
{ 0 ,1/ ( p - 1) ,2/ ( p - 1) , . . . ,1} 中随机取值 ;
3) 设 P 3 为 k ×k 的随机混淆矩阵 (random permutation
matrix) ,矩阵中的每一列及每一行都只有一个值为 1 ,其余
值都为 0 .
显然 ,D 3 、X3 和 P 3 中随机取值都是相互独立的. 令B 3
= (J m ,1 X3 + ( △/ 2) [ (2B - J m , k) D 3 + J m , k ]) P 3 为矩阵B的
随机化矩阵. 和 △B 相似 ,B 3 中的每相邻两行中也只有一个
参数的取值不同 ,每相邻 两行作为模型参数输入 ,根据方程
(2)可获得模型每一个参数的灵敏度. 但是 ,由于 Morris 法的
随机性 ,很容易在一次随机取样及随机化过程中出现误差 ,
所以可以进行多次重复 ,取多次平均值表征所有参数的灵敏
度 ,确定其大小排序 ;而其标准差可用来表征参数之间相互
作用的程度 ,如果标准差小 ,说明该参数与其它参数之间的
相互作用程度小 ;如果标准差大 ,则说明该参数与其它参数
的相互作用程度大. Morris 法在确定模型各参数灵敏度大小
排序时简单而有效 ,可用来“冻结”那些灵敏度很小的模型参
数 ,对剩下的参数再做定量的全局灵敏度分析. Crosetto 和
Tarantola[7 ]对一个基于 GIS 的洪水预测模型进行灵敏度分
析时先用 Morris 法筛选出对模型结果影响最大的 5 个参数 ,
再用定量全局灵敏度分析的方法量化 (傅里叶幅度灵敏度检
验扩展法) 5 个参数中的每一个参数的不确定性对模型结果
的影响.
313 傅里叶幅度灵敏度检验法
傅里叶幅度灵敏度检验法 ( Fourier Amplitude Sensitivity
Test)是在 20 世纪 70 年代引入模型的灵敏度分析的[10~12 ] ,
是那个时代甚至是当今灵敏度分析中最好的方法之一. 该方
法的核心是用一合适的搜索曲线在参数的多维空间内搜索 ,
从而把多维的积分转化为一维的积分. 该方法对模型中每一
个输入参数都引入一个具有共同独立参数 ( s) 的函数 ,并给
每一个输入参数定义一个整数频率 ,最后使模型成为独立参
数的周期函数. 对模型的输出结果进行傅里叶分析 ,产生每
一个频率的傅里叶幅度. 用傅里叶幅度的大小来指示每一个
参数的灵敏度 ,幅度越大 ,说明模型对该参数越敏感 ;反之亦
然 ,该幅度也被称为重要性测度 [21 ,22 ,24 ,41 ] .
8501 应 用 生 态 学 报 15 卷
对于生态模型 y = f ( x1 , x2 , ⋯, xk) ,其中 , k 是模型中
的参数个数 , FAST 定义单个参数的灵敏度为由该参数的微
小变化所引起的在整个输入参数空间内的模型结果的平均
变化 :5 y5 x i =∫x 5 y ( X)5 x i P( X) dX (3)
其中 P( X) 为 X = ( x1 , x2 , ⋯, xk) 的联合分布概率. 设
x i = gi (sin ( w is) ) , i = 1 ,2 , ⋯, k (4)
w i 是一个人为设定的频率 , s 是所有参数的独立参数 , gi 称
为搜索函数. 由 (4) ,方程 (3) 可转化为 ,5 y5 x i =∫s 5 y ( x1 ( s) , x2 ( s) , ⋯, xk ( s) )5 x i ( s) P( x1 ( s) , x2 ( s) ,
⋯, xk ( s) ) ds (5)
k 维积分转化为一维积分. 为了保证当 s 变化时 ,由方程 (4)
获得的 x i 与 x i 的概率分布一致 ,必须满足如下方程组 :
π 1 - u2 f i ( gi ( u) )
d gi ( u)
d u = 1 , gi (0) = 0 (6)
如果 w i 非线性相关 ,而函数 gi 又满足条件 (6) ,那么 (5) 和
(3)的结果相等. 所谓非线性相关 ,需满足条件6k
i = 1
a i w i ≠0 (7)
其中 ,a i 为整数. 如果使用非线性相关的频率 ,那么方程 (5)
的积分计算需要 s 在无限的空间内取值 ,这在计算上是不可
行的. 如果采用整数的频率 ,那么方程 (5)的积分计算只要求
s 在一个有限的区间 ( [ - π,π])内取值. 这样 ,模型就可以表
达为
f ( s + 2π) = f ( s) (8)
这是一个关于 s 的周期性函数 ,可展开为傅里叶级数
f ( s) = 6k
i = 1
A i sin ( w is) (9)
其中 A i 由如下方程获得
A i =
1
π∫
π2πf ( s) sin ( w is) ds (10)
用 A i 表示第 i 个参数的灵敏度. 对 s 在区间[ - π,π] 内等间
隔取样 ,把取样获得的每一个参数输入模型 ,多次运行模型 ,
由如下方程可近似获得 A i ,
A i =
α
Ns 6Nsk = 1 f ( sk) sin ( w isk) (11)
其中 N s 为取样数. Cukier 等[10 ]建议
N s = 2M w max + 1 (12)
对于搜索函数 gi ,如果 x i 服从[0 ,1 ]的均匀分布 ,Cukier 等[10 ]
建议
x i =
x ie
v isin ( w is) (13)
其中 ,
x i 为中心值 ,
v 为 x i 区间的端点值 , s 在[ - π/ 2 ,π/ 2 ]
内变化. Koda 等[26 ]建议
x i =
x i (1 +
v i sin ( w is) ) (14)
Saltelli 等[43 ]指出 ,上述两个搜索函数并不能充分在参数空
间内搜索 ,并提出新的搜索函数
x i =
1
2 +
1
πarcsin (sin ( w is) ) (15)
上述的搜索函数都只是针对参数是 [ 0 ,1 ]均匀分布 ,对于分
布函数为 F i ( x) 的参数 ,Lu 等[27 ]指出 ,其搜索函数应为
x i = F - 1i ( c + 1πarcsin (sin ( w is) ) ) (16)
其中 , F - 1i ( x) 为参数 i 的分布函数的反函数 , c 为常数.
上述 3 种方法是定性全局灵敏度分析中比较常用的方
法 ,还有很多其它的方法 , 如析因实验 (factorial experi2
ment) [20 ]和 Plackett2Burman 法[2 ,9 ] .
4 定量全局灵敏度分析
定量全局灵敏度分析定量地确定各模型参数对于模型
结果中误差的贡献率 . 其主要方法有 Sobol’法[47 ]和傅里叶
幅度灵敏度检验扩展法 ( Extended Fourier Amplitude Sensi2
tivity Test) [43 ] . 这两种方法都是基于方差的方法 ,认为模型
结果的方差可完全反映模型结果的不确定性. 它们不单单计
算参数对模型结果的单独影响 ,还考虑参数之间的相互作用
对模型结果的影响. 在做定量全局灵敏度分析时 ,可以先做
定性的全局灵敏度分析 ,从而过滤一些对模型结果影响不大
的参数[3 ] .
411 Sobol’法
Sobol’法的核心是把模型分解为单个参数及参数之间
相互组合的函数. 假设模型为 Y = f ( x) , ( x = x1 , x2 , ⋯,
xk) , x i 服从[0 ,1 ] 均匀分布 , f 2 ( x) 可积. 模型可分解为
f ( x) = f 0 + 6k
i = 1
f i ( x i) + 6
i < j
f ij ( x i , x j) + ⋯ +
f 1 ,2 , ⋯, k ( x1 , x2 , ⋯, xk) (17)
方程右边共有 2 k 项. 但是这种分解并不是唯一的. 但如果满
足下列条件
∫
1
0
f i ( x i) dx i = 0 , Π x i , i = 1 ,2 , ⋯, k
∫
1
0
1
0
f ij ( x i , x j) dx i dx j = 0 , Π x i , x j , i < j (18)
∫
Ω
f 12 . . . k ( x1 , x2 , ⋯, xk) dx1 dx2 ⋯dxk = 0
则方程 (17) 具有唯一的分解形式. 分解的各项满足
∫Ωkf i1 , ⋯, isf j1 , ⋯, jl dx = 0 , ( i1 , ⋯, is) ≠( j1 , ⋯, jl) , k = s + l
(19)
记
∫Ωf ( x) dx = f 0 (20)
对除 x i 外的所有参数积分可获得 f i ( x i)
∫f ( x) 7j ≠i dx j = f 0 + f i ( x i) (21)
对除 x i , x j 外的所参数积分可获得 f ij ( x i , x j)
∫f ( x) 7l ≠i , jdx l = f 0 + f i ( x i) + f j ( x j) + f ij ( x i , x j)
(22)
依此类推 , 可获得 f 1 ,2 , ⋯k ( x1 , x2 , ⋯, xk) . Sobol’用总的方
差
95016 期 徐崇刚等 :生态模型的灵敏度分析
V =∫f 2 ( x) dx - f 20 (23)
来表示所有参数对模型结果的影响程度 ,用偏方差
V i =∫f 2i dx i (24)
来表示单个参数对模型结果的影响程度 ,用偏方差
V i1 , i2 , ⋯, is =∫f 2i1 , ⋯, is dx i1 dx i2 ⋯dx is (25)
来表示参数之间的作用对模型结果的影响程度. 对方程 ( 17)
两边平方再积分可得
V = 6k
i = 1
V i + 6
i < j
V ij + ⋯+ V 1 ,2 , ⋯k (26)
对上式各项归一化 ,并令
S i1 , ⋯, is =
V i1 , ⋯, is
V
(27)
可获得模型单个参数及参数之间相互作用的灵敏度 S. 于
是 ,方程 (26) 可改写为
1 = 6k
i = 1
S i + 6
i < j
S ij + 6
i < j <1
Sijl + ⋯+ S1 ,2 , ⋯, k (28)
对于 Si , 称之为一次灵敏度 ; S ij 为二次灵敏度 , 依此类推 ,
S12 . . . k 为 k 次灵敏度. 引入参数 i 的总灵敏度 S Ti
S Ti = 6 S ( i) (29)
S( i) 指所有包含参数 i 的灵敏度. 因为方程右边一共有 2 k - 1
项 ,如果 k 很大则很难实现. 因此 ,考虑更一般的情况 ,如果
我们需要检验模型对一组参数的灵敏度 S Ti ,那么 ,可以先把
所有参数分成两组 , u 和 v . 这样
x = ( x1 , x2 , ⋯, xk) = ( u , v)
f ( x) = f 0 + f 1 ( u) + f 2 ( v) + f 12 ( u , v)
V =∫f 2 ( x) dx - f 20 (30)
V 1 =∫f 21 ( u) du
V 2 =∫f 22 ( v) dv
V 12 = V - V 1 - V 2
根据方程 (21) , f 1 ( u) =∫f ( u , v) dv - f 0 ,所以
V 1 = (∫f ( u , v) dv) 2 - f 20
=∫f ( u , v) dv∫f ( u , v′) dv′- f 20 (31)
采用蒙特卡罗法
V 1 + f 20 ≈ 1N 6Nj = 1 f ( uj , vj) f ( uj , v′j)
f 0 ≈ 1N 6Nj = 1 f ( uj , vj) (32)
V + f 20 ≈ 1N 6Nj = 1 f 2 ( uj , vj)
V 2 + f 20 ≈ 1N 6Nj = 1 f ( uj , vj) f ( u′j , v j)
其中 uj , v j 和 u′j , v′j 分别为两次独立蒙特卡罗取样的一个
样本 , N 为每次蒙特卡罗取样的样本数量. 假设 k = 4 ,而要
计算 S 3 = V 3/ V ,那么只需要把模型的参数分成两组 u =
( x1 , x2 , x4) , v = x3 对模型的各参数分两次取样 (每次取样
的样本数量为 N) ,由方程 (32) 可知
V 3 + f 20 ≈ 1N 6Nj = 1 f ( x j1 , x j2 , x j3 , x j4) f ( x′j1 , x′j2 , x j3 ,
x′j4) (33)
其中 , ( x j1 , x j2 , x j3 , x j4) 及 ( x′j1 , x′j2 , x′j3 , x′j4) 分别为第一
和第二次取样时的样本. 而 S T3 可由如下方程组获得
S T3 = 1 - V - 3/ V
V
- 3 + f 20 ≈ 1N 6Nj = 1 f ( x j1 , x j2 , x j3 , x j4) f ( x j1 , x j2 , x′j3 ,
x j4) (34)
显然 ,Sobol’通过对模型结果的方差的分解 ,定量地获
得每一个参数的一次及高次灵敏度. 如果计算各参数的总敏
度 S Ti ,通过归一化就可获得每一参数的相对贡献率. Sobol’
法已开始应用于生态模型的灵敏度分析. Pastres 等[37 ] 用
Sobol’法对一个三维浅水富营养化模型进行定量全局灵敏
度分析 ,结果表明 ,含氮量是控制整个水生生态系统第一生
产力的主要因素 ,但是深水藻类的初始状态是导致“水华”产
生的更重要的因素 ; Salvador 等[44 ]用 Sobol’法对 Rothermel
火传播模型在地中海灌木林区应用的 10 个参数进行全局灵
敏度分析 ,结果表明 ,低热含量、矿物含量及细粒密度等 3 个
参数对模型输出结果 (火线上某一点的传播速度和火烧产生
的能量)的影响很小 ,而其它 7 个参数的影响都较显著. 此
外 ,Sobol’法中对参数分组的思想也在生态模型的局部灵敏
度分析中得到了应用. KÊhler 和 Wirtz[25 ]对一个水生生态系
统模型 (AQ EM) 进行灵敏度分析时 ,先把模型参数分成几
组 ,把一组参数看成一个参数. 组内参数的变化按如下规则 :
1)每一个参数都增加 (或减小) 50 % ;2) 每一个参数增加 (或
减小)两个偏差. 这样就大大减小了局部灵敏度分析的工作
量.
412 傅里叶幅度灵敏度检验扩展法
Saltelli 等结合 Sobol’法和傅里叶幅度灵敏度检验法的
优点 ,提出了傅里叶幅度灵敏度检验扩展法 ( Eextended
Fourier Amplitude Sensitivity Test) [43 ] . 该方法由傅里叶转换
获得傅里叶级数的频谱 ,通过该频谱曲线分别得到由每一个
参数及参数的相互作用所引起模型结果的方差.
根据合适的搜索函数 ,模型 y = f ( x1 , x2 , ⋯, xk) 可转化
为 y = f ( s) . 通过傅里叶变换
y = f ( s) = 6+ ∞
j = - ∞
{ A jcos js + B j sin js} (35)
其中
A j =
1
2π∫
π
- πf ( s) cos jsds ,
B j =
1
2π∫
π
-
πf ( s) sin jsds
(36)
j ∈Z = { - ∞, ⋯, - 1 ,0 , 1 , ⋯, + ∞} . 傅里叶级数的频谱曲
线定义为Λj = A2j + B2j ,其中 , j ∈Z ,A - j = A j ,B - j = B j ,Λ- j
=Λj . 由参数 x i 不确定性所引起的模型结果的方差
V i = 6
p ∈Z0
Λpw i = 2 6+ ∞
p = 1
Λpw i (37)
其中 Z0 = Z - { 0} . 总的方差为
0601 应 用 生 态 学 报 15 卷
V = 6
j ∈Z0
Λj = 2 6+ ∞j = 1Λj (38)
对 s 在区间[ - π,π]内等间隔取样 ,把取样获得的每一个参
数输入模型 ,多次运行模型 ,由如下方程可近似获得 A j 和
B j
A j =
1
Ns 6Nsk = 1 f ( sk) cos( jsk)
B j =
1
Ns 6Nsk = 1 f ( sk) cos( jsk) , j ∈
Z (39)
其中 , Ns 为取样数 ,
Z = { - Ns - 12 , ⋯, - 1 , 0 , 1 , ⋯, Ns - 12 }< Z , Ns 满足方程 (12) . 由 A i 和 B j 及参数 x i 所对应的频率
w i ,通过方程 (37) 、(38) 可获得每一个参数所引起的方差 V i
及模型结果的总方差 V . 通过方程 (27) ,可分别获得每一个
参数的灵敏度. 要计算参数 x i 的总灵敏度 (见方程 (29) ) ,可
以先给 x i 设定一个频率 w i ,而为剩余的所有其它参数设定
一个不同的频率 w′i . 通过计算 w′i 在 pw′i 上的所有频谱
值 ,就可得到偏方差 V
- i ,它包含除 x i 外的所有参数及其相
互关系的影响. 所以参数 x i 的总灵敏度 S Ti = ( V - V - i ) /
V . 对各参数的总灵敏度归一化 ,就可获得参数 x i 的不确定
性对模型结果总的不确定性的贡献率.
EFAST 的优点在于其计算量要比 Sobol’法小得多.
Sobol’法需要两次取样才能得到某个参数的一次灵敏度和
总灵敏度 ,而 EFAST 只需要一次取样就可同时获得某个参
数的一次灵敏度和总灵敏度. EFAST 和 Sobol’法都要求模
型的各输入参数是不相关的. 目前 EFAST 已在模型的灵敏
度分析上获得广泛应用 [6~8 ,48 ] ,但在生态模型的灵敏度分析
中的应用还未见报道.
5 展 望
生态预测是生态学新近出现的一个新领域 [5 ] . 生态模型
在其中扮演着非常重要的角色. 随着生态模型的发展 ,生态
模型的灵敏度分析也愈发显得重要. 但是 ,目前使用的生态
模型灵敏度分析方法比较单一 ,以局部灵敏度分析为主 ,开
始大量出现定性的全局灵敏度分析 ,但是主要是采用基于拉
丁超几何体采样的多元回归法. 对于另一种常见的定性全局
灵敏度分析法 (傅里叶幅度灵敏度检验法)则未见报道. 这可
能与傅里叶幅度灵敏度检验法比多元回归法较复杂有关. 定
量全局灵敏度分析是当今模型灵敏度分析的前沿和热点 ,但
在生态模型中的应用才刚刚起步. 对于国内生态模型的灵敏
度分析 ,大多局限于局部灵敏度分析 ,未见全局灵敏度分析
的报道. 结合国内外生态模型和灵敏度分析方法的发展现
状 ,提出了今后生态模型灵敏度分析的热点和难点.
511 联合灵敏度
定量全局灵敏度分析是当今模型灵敏度分析的前沿和
热点. 但现有的两种定量全局灵敏度分析方法 (傅里叶幅度
灵敏度检验扩展法和 Sobol’法)都只能定量获得某一个参数
的单独灵敏度及总灵敏度 ,而不能确定两个或多个参数联合
灵敏度. Crosetto 和 Tarantola[7 ]指出 ,通过比较 S i 和 S Ti可以
知道单个参数对模型结果的影响和该参数与其参数联合作
用对模型结果影响的差别. 但这并不能得到具体参数之间的
联合灵敏度. 在生态学模型中 ,定量两个参数之间的相互作
用对模型结果的影响有很重要的生态学意义. 比如 ,火烧和
采伐都会单独影响森林景观动态 ,但它们之间的相互作用也
会在很大程度上影响景观动态. 采伐会减少可燃物的积累从
而会降低火烧强度 ;而火烧又会减少木材量 ,从而增加采伐
面积 (假设采伐量是一定的) . 因此 ,对景观动态模型进行灵
敏度分析时 ,除了单独考虑火烧和采伐两个因子 ,还应该考
虑两个因子间的相互作用.
512 空间直观景观模型的灵敏度分析
现有的灵敏度分析法均建立在传统模型或者说非空间
模型之上. 随着空间直观景观模型的大量涌现 ,怎样把传统
的灵敏度分析法应用到空间直观景观模型上是一亟待解决
的问题. 大部分空间直观景观模型是建立在栅格数据的基础
上 ,把景观概念化为一由相同大小的象元或样地组成的格
网[32 ] .由于一般研究区都是由成千上万个样地组成 ,灵敏度
分析不可能计算每一样地内的不确定性对模型结果的影响.
Crosetto 等[6 ]在研究二维空间参数的不确定性对模型结果的
影响时 ,通过加入或不加入每个象元的随机误差所引起模型
结果的变化程度来定量化该空间参数的灵敏度. 这是一种很
好的解决方法 ,但其难点在于怎样在每次取样时产生该空间
参数的随机误差 ,即怎样建立空间参数的误差模型.
513 模型灵敏度的共性
几乎所有的模型都是基于具体实例的 ,从一个模型得到
的结果不能外推到其它的情况. 也许最直接的解决方法是设
计出相对通用的生态模型. 模型的共性可以通过比较几个模
型来总结 ,如通过对一系列的模型进行灵敏度分析可以发现
模型对一些参数的灵敏度较差 ,那么这个低灵敏度可以认为
是这组模型的共性 [14 ] . 但是目前对此还没有深入的研究.
参考文献
1 Battaglia M ,Sands P. 1998. Application of sensitivity analysis to a
model of Eucalypt us globulus plantation productivity. Ecol Model ,
111 :237~259
2 Beres DL ,Hawkins DM. 2001. Plackett2Burman technique for sen2
sitivity analysis of many2parametered models. Ecol Model ,141 :171
~183
3 Campolongo F ,Saltelli A. 1997. Sensitivity analysis of an environ2
mental model :An application of different analysis methods. Reliab
Engin Syst S af ety ,57 :49~69
4 Carlson DH , Thruow TL . 1996. Comprehensive evaluation of im2
proved SPUR model (SPUR~91) . Ecol Model ,85 :229~240
5 Clark J S , Carpenter SR ,Barber M , et al . 2001. Ecological fore2
casts :An emerging imperative. Science ,27 :657~660
6 Crosetto M ,Ruiz JAM ,Crippa B. 2001. Uncertainty propagation in
models driven by remotely sensed data. Rem Sens Envi ron ,76 :373
~385
7 Crosetto M , Tarantola S. 2001. Uncertainty and sensitivity analy2
sis : Tools for GIS2based model implementation. IN T J Geogr Iin2
f orm Sci ,15 (5) :415~437
8 Crosetto M , Tarantola S ,Saltelli A. 2000. Sensitivity and uncertain2
ty analysis in spatial modeling based on GIS. A gric Ecosyst Envi2
ron ,81 :71~79
9 Cryer SA , Havens , PL . 1999. Regional sensitivity analysis using a
16016 期 徐崇刚等 :生态模型的灵敏度分析
fractional factorial method for the USDA model GL EAMS. Envi2
ron Model Sof tw are ,14 :613~624
10 Cukier RI ,Fortuin CM ,Shuler KE , et al . 1973. Study of the sensi2
tivity of coupled reaction systems to uncertainties in rate coefficients
I. Theory. J Chem Phys ,59 :3873~3878
11 Cukier RI ,Levine HB ,Shuler KE. 1978. Nonlinear sensitivity anal2
ysis of multiparameter model systems. J Com put Phys ,26 :1~42
12 Cukier RI ,Schaibly J H ,Shuler KE. 1975. Study of the sensitivity of
coupled reaction systems to uncertainties in rate coefficients Ⅲ.
Analysis of the approximations. J Chem Phys ,63 :1140~1149
13 Downing DJ , Gardner RH , Hoffman FO. 1985. An examination of
response2surface methodologies for uncertainty analysis in assess2
ment models. Technomet rics ,27 :151~163
14 Dunning JB ,Stewart DJ ,Danielson BJ , et al . 1995. Spatially explic2
it population models :Current forms and future uses. Ecol A ppl ,5 :3
~11
15 Fennel K ,Losch M , Schroter J , et al . 2001. Testing a marine e2
cosystem model : Sensitivity analysis and parameter optimization. J
M ar Syst ,28 :45~63
16 Gardner RH ,O’Neill RV ,Mankin JB , et al . 1981. A comparison of
sensitivity analysis and error analysis based on a stream ecosystem
model. Ecol Model ,12 :173~190
17 Gentil S ,Blake G. 1981. Validation of complex ecosystem models.
Ecol Model ,14 :21~38
18 Halvorsen E ,Pedersen OP ,Slagstad D , et al . 2001. Microzooplank2
ton and mesozooplankton in an upwelling filament off Galicia :Mod2
elling and sensitivity analysis of the linkages and their impact on the
carbon dynamics. Prog Oceanogr ,51 :499~513
19 Helton J C. 1993. Uncertainty and sensitivity analysis techniques for
use in performance assessment for radioactive waste disposal. Reliab
Engin Syst S af ety ,42 :327~367
20 Henderson2Sellers B ,Henderson2sellers A. 1996. Sensitivity evalua2
tion of environmental models using fractional factorial experimenta2
tion. Ecol Model ,86 :291~295
21 Homma T ,Saltelli A. 1996. Importance measures in Global sensitiv2
ity analysis of Nonlinear models. Reliab Engin Syst S af ety ,52 :1~
17
22 Iman RL , Hora SC. 1990. A robust measure of uncertainty impor2
tance for use in fault tree system analysis. Risk A na ,10 :401~406
23 Iman RL . 1999. Latin hypercube sampling. Encyclopedia of Statisti2
cal Science. Vol. 3. New York :Wiley. 408~411
24 Ishigami T , Homma T. 1990. An importance quantification tech2
nique in uncertainty analysis for computer models. Proceedings of
the ISUMA’90 , First International Symposium on Uncertainty
Modelling and Analysis ,University of Maryland ,USA. New York :
IEEE Computer Society. 398~403
25 KÊhler P ,Wirtz KW. 2002. Linear understanding of a huge aquatic
ecosystem model using a group2collecting sensitivity analysis. Envi2
ron Model Sof tw are ,17 :613~625
26 Koda M ,McRae GJ ,Seinfeld J H. 1979. Automatic sensitivity analy2
sis of kinetic mechanisms. Intern J Chem Kinet ,11 :427~444
27 Lu YC , Mohanty S. 2001. Sensitivity analysis of a complex ,pro2
posed geologic waste disposal system using the Fourier Amplitude
sensitivity test method. Reliab Engin Syst S af ety ,72 :275~291
28 Mahamah DS. 1988. Simplified sensitivity analysis applied to a nu2
trient2biomass model. Ecol Model ,42 :103~109
29 Majkowski J , Tidgeway J M ,Miller DR. 1981. Multiplicative sensi2
tivity analysis and its role in development of simulation models. Ecol
Model ,12 :191~208
30 McCarthy MA ,Burgman MA , Ferson S. 1995. Sensitivity analysis
for models of population viability. Biol Conserv ,73 :93~100
31 Mc Kay MD , Beckman RJ , Conover WJ . 1979. A comparison of
three methods fro selecting values of input variables in the analysis
of output from a computer code. Technomet rics ,21 :239~245
32 Mladenoff DJ , He HS. 1999. Design and behavior of LANDIS ,an
object2oriented model of forest landscape disturbance and succes2
sion. In :Mladenoff DJ ,Baker WL ,eds. Spatial Modeling of Forest
Landscape Change :Approaches and Applications. Cambridge :Cam2
bridge University Press. 125~162
33 Morris MD. 1991. Factorial sampling plans for preliminary compu2
tational experiments. Technomet rics ,33 :161~174
34 Nathan R ,Safriel UN ,Noy2Meir I. 2001. Field validation and sensi2
tivity analysis of a mechanistic model for tree seed dispersal by
wind. Ecology ,82 (2) :374~388
35 Omlin M ,Brun R ,Reichert P. 2001. Biogeochemical model of Lake
Zurich : Sensitivity , identifiability and uncertainty analysis. Ecol
Model ,141 :105~123
36 Pacala ,SW ,Canham CC ,Saponara J , et al . 1996. Forest models de2
fined by field measurement : Estimation ,error analysis and dynam2
ics. Ecol Monogr ,66 :1~43
37 Pastres R ,Chan K ,Solidoro C , et al . 1999. Global sensitivity analy2
sis of a shallow water 32D eutrophication model. Com put Phys
Com m un ,117 :62~74
38 Pastres R ,Franco D ,Pecenik G , et al . 1997. Local sensitivity analy2
sis of a distributed parameters water quality model. Reliab Engin
Syst S af ety ,57 :21~30
39 Ratto M , Tarantola S ,Saltelli A. 2001. Sensitivity analysis in model
calibration : GSA2GLU E approach. Com put Phys Com m un , 136 :
212~224
40 Recknagel F. 1984. A comprehensive sensitivity analysis for an eco2
logical simulation model. Ecol Model ,26 :77~96
41 Saltelli A , Anderes TH , Homma T. 1993. Sensitivity analysis of
model output :An investigation of New Techniques. Com put S tatist
Data A nal ,15 :211~238
42 Saltelli A ,Chan K ,Scott M. 2000. eds. Sensitivity Analysis , Proba2
bility and Statistics Series. New York :John Wiley & Sons.
43 Saltelli A , Tarantola S ,Chan KPS. 1999. A quantitative model2inde2
pendent method for global sensitivity analysis of model output .
Technomet rics ,41 :39~56
44 Salvador R , Pinnol J , Tarantola S , et al . 2001. Global sensitivity
analysis and scale effects of a fire propagation model used over
Mediterranean shrublands. Ecol Model ,136 :175~189
45 Sharov AA. 1995. Non2linear sensitivity analysis as a tool for ex2
plaining population dynamics in a life2system conceptual frame2
work. Com put Elect ron A gric ,13 :103~113
46 Shirley MDF , Rushton PS , Yong AG , et al . 2001. Simulating the
long2term dynamics of slug populations : A process2based modeling
approach for pest control. J A ppl Ecol ,38 :401~411
47 Sobol’IM. 1993. Sensitivity Estimates for nonlinear mathematical
models. M ath Model Com put Ex p ,1 :407~414
48 Tarantola S , Giglioli N ,Jesinghaus J , et al . 2002. Can global sensi2
tivity analysis steer the implementation of models for environmental
assessments and decision2making ? S tochastic Envi ron Res Risk
Asses ,16 :63~76
49 Tomovic′R. 1963. Sensitivity Analysis of Dynamic Systems. New
York :Mc Graw2Hill. 142
50 Waide JB ,Webster J R. 1976. Engineering systems analysis. In : Pat2
ten BC ,eds. Systems Analysis and Simulation in Ecology. 4. New
York :Academic Press. 329~371
51 Zhang H (张 宏) , Fang Z2L (樊自立) . 2000. Studies on a NPP
model of salinized meadow in the north of Tarim basin. Acta Phy2
toecol S in (植物生态学报) ,24 (1) :13~17 (in Chinese)
52 Zhu Z2R(祝增荣) , Hang C2W (黄次伟) ,Chen J2A (程家安) , et
al . 1994. Simulation analysis of white2backed flanthopper popula2
tion dynamics on first season rive in Zhejiang Province ,China. Acta
Ecol S in (生态学报) ,14 (2) :188~195 (in Chinese)
作者简介 徐崇刚 ,男 ,1979 年生 ,博士生 ,主要研究方向为
空间直观景观模型的不确定性和灵敏度分析 ,发表论文 5
篇. E2mail :xuchongang @yahoo. com
2601 应 用 生 态 学 报 15 卷