全 文 ：林木生长的多维时间序列分析 3
吴承祯 3 3 洪 伟 (福建林学院资源与环境学系 ,南平 353001)
【摘要】 利用多维时间序列分析方法 ,以影响杉木直径生长的五大主导气象要素作控制因子 ,建立杉木直径生
长的 CAR 模型 ,从而对杉木直径生长提前一年作预测 ,回验结果表明模型准确率很高 ,为林木生长预测预报提
供了一种新方法.
关键词 CAR 模型 生长过程
Multidimensional time series analysis on tree growth. Wu Chengzhen and Hong Wei ( Depart ment of Resources and
Envi ronment , Fujian Forest ry College , N anping 353001) . 2Chin. J . A ppl . Ecol . ,1999 ,10 (4) :395～398.
By the analysis of multidimensional time series , five important climatic factors were used as control factors in CAR
model to forecast the diameter growth of Chinese fir. This model could predict the diameter growth of Chinese fir one
year in advance , and the accuracy of forecasting results was quite high. This method may be regarded as a new way for
predicting tree growth.
Key words CAR model , Growth process.
3 福建省自然科学基金资助项目 ( F991) .
3 3 通讯联系人.
1997 - 06 - 23 收稿 ,1998 - 03 - 30 接受.
1 引 言
时间序列分析是时间数列预测和回归预测方法的
综合形式 ,具有较高的预测准确性 ,并适用于多种复杂
条件下的预测. 正因为此 ,时间序列分析已广泛应用于
农、林业生产的生长、病虫害、土壤性质变化等领
域[1 ,4～11 ,14 ,18 ] . 前人[4 ,7 ,8 ,14 ]在应用时间序列方法建立
树木生长量模型时 ,都采用一维时间序列 ,但林木生长
受环境等生态因素的影响 ,其处于一个开放的系统 ,因
此必须在一维时间序列模型中引入控制项 ,即用多维
时间序列来研究林木生长预测问题. 多维时间序列分
析是系统工程学上为解决动态系统的控制问题而提出
的一种建模方法 ,其主要特征是在考虑因变量时序变
动的基础上 ,融入了自变量的控制作用.
林木生长的数学模拟可用于预测生长量 ,也是分
析林木生长动态的有力手段. 长期以来 ,人们曾提出包
括一维时间序列模型在内的各种各样的数学模型对林
木生长过程进行定量描述. 本文将多维时间序列分析
建模方法应用于林木生长预测和动态分析 ,由于引入
影响林木生长的气象生态因子为控制因子 ,该方法将
为解决林木生长预测问题提供一个新的最优化研究方
法 ,也将丰富树木气候学研究内容.
2 材料与方法
2. 1 资料的收集
根据全国杉木产区区划 ,福建尤溪县属于一般产区 ,其杉
木地理分布范围位于中心产区边缘. Cox 认为 [17 ] ,靠近树种地
理分布或垂直分布边缘地区 ,更能提供对各种气候因子年敏感
的生长材料. 为此 ,选定尤溪县 120 株杉木优势木解析木材料 ,
以年为单位测定其年轮宽度 ,其测定精度要求 0. 1cm.
决定杉木年轮宽度的因子 ,除气候条件外 ,还受树龄因子
的影响. 采用年轮指数 ( Y) 作为因变量以消除树龄的影响 ,其
计算公式为[2 ] :
Y = 该年实际年轮宽度/ 树龄因子生长的年轮宽度
式中 ,分母值是建立该年实际年轮宽度 ( y) 与树龄 ( x ) 的指数
曲线方程 y = abx ,尔后将各树龄 ( x)值分别代入方程求出的回
归计算值.
气候因子的选择 ,根据杉木生物学特性 ,选择全年日平均
气温在 18～27 ℃的日数、5 月和 6 月平均最小相对湿度、15cm
年平均地温等因子.
2. 2 多维时间序列模型
众所周知 ,一阶自回归 AR 模型形式为 :
X t =θ1 X t - 1 + ⋯+θnX t - n +εt (1)
此模型即为前人研究时间序列常用的时间序列模
型[1 ,7～11 ,14 ,17 ] .
n 阶自回归滑动平均 ARMA 模型形式为 :
X t - α1 X t - 1 - ⋯- αnX t - n =εt - β1εt - 1 - ⋯- βnεt - n (2)
即文献[10 ]所采用方法.
如果在式 (2)中 ,再加入受控制项 ,则有如下形式 :
X t - α1 X t - 1 - ⋯- αnX t - n - β0 U t - β1 U t - 1 - ⋯- βnU t - n
=εt - γ1εt - 1 - ⋯- γnεt - n (3)
式 (3)称为 CARMA 模型 ,即为带控制项的自回归滑动平
均模型 ,也就是经典的多维时间序列分析模型.
由于 CARMA 模型十分复杂 ,邓自立[3 ] 提出了特殊的
应 用 生 态 学 报 1999 年 8 月 第 10 卷 第 4 期
CHIN ESE JOURNAL OF APPL IED ECOLO GY ,Aug. 1999 ,10 (4)∶395～398
CARMA 模型 ,即用 CAR (带受控项的自回归) 模型 ,而对动态
系统实行统一建模的新观点. 两者之间的关系是 :任何 CARMA
模型均可以充分高阶的 CAR 模型逼近到任意精度 ,因而实际
上可用 CAR 模型对动态系统统一建模. 通常用极大似然 (ML)
方法辩识 CARMA 模型 ,可得到模型参数的一致估计 ,但算法
复杂 ,计算量大. 虽然用递推增广最小二乘法 ( REIS) 建立
CARMA 模型是简单的 ,但参数估计可能有偏. 用 CAR 模型对
系统建模的优点是 :用普通递推最小二乘法 (RLS) 可得到模型
参数的一致估计 ,而且计算简单 ,便于在微机上实现 ,克服了用
CARMA 模型对动态系统建模在计算机上的缺陷.
CAR 模型自动辩识机理是从低阶开始对系统建模 ,然后逐
次增加模型阶数 ,并用 F 检验对这些模型进行自动筛选. 它由
3 部分组成 :
1)递推最小二乘法估计器
假定动态系统可用如下 CAR 模型描写 :
Y t =α1 Y t - 1 + ⋯+αn Y t - n + b0 U t + ⋯+ bnU t - n +εt (4)
式中 , U 为控制项 ,利用数理统计中的递推最小二乘法 [16 ] ,可
估计出上式中的各个变量参数及残差εt ,控制项可以为多个.
2)最高阶 n 的 F 检验判决器
辩识机所采用的 CAR 模型定阶法是基于已知的 N 对观
测数据 ( U t , Y t , t = 1 ,2 , ⋯, N ) ,由低阶到高阶递推地对系统拟
合 CAR 模型并依次对相邻的两个 CAR 模型用 F 检验法判别
模型阶的合适性. 对于 3 因素控制的 CAR 模型 ,在置信度 1 - a
= 5 %时 ,通常取 Fa = 3 ,即当相邻的 F 值小于 Fa 时 ,CAR ( n)
模型是合适的. 假定真实值模型的阶是 n ,可由低阶开始递推
地建立 CAR (1) ,CAR (2) , ⋯,CAR ( n - 1) ,CAR ( n) , CAR ( n
+ 1) ,且相邻的每两个模型间用上述 F 检验法检验 ,直至找到
合适的模型 CAR( n + 1)为止 ,一共只需建立 n + 1 个模型.
3)真实阶及时滞的 F 检验判决器
虽然得到了合适的 CAR( n) 模型 ,但其中某些参数可能接
近于零 ,因而可用 F 检验法把这些变量剔除掉 ,如果式 (4) 中
b0 = b1 = b2 = ⋯bk - 1 = 0 ,则系统的时滞为 k .
从所得到的 CAR( n)中剔除掉接近于零的参数后 ,重新用
递推最小二乘法建立参数的 CAR( n)模型. 用 F 检验法对原来
拟合的 CAR( n)模型和这个新建立的较少参数的 CAR( n)模型
进行检验 ,如不显著 ,则较少参数的 CAR ( n) 模型是真实的 ;如
显著 ,则原来的 CAR( n)模型为真实的.
2. 3 突变预测模型
含有突变因素的 GM (1 ,1) 模型反映了环境突变因素的影
响和作用 ,其灰微分影子方程为
dx
(1)
dt + ax (1) = u
式中 ,a 和 u 分别为环境因素突变后 ,系统行为新的发展系数与
新的灰作用量 ,含有突变因素的 GM (1 ,1)模型的解为 :
X^ (1)( k + 1) = ( X (0)(1) - ua ) ·e - ak + ua (5)
X^ (0)( k + 1) = ( X (0)(1) - ua ) ·e - ak (1 - ea) (6)
模型 (5)由如下序列建模生成 ,即〔x ( k) │k = 1 , 2 , ⋯, 2 +
q〕
=〔x (1) , x (2) , ⋯, x ( q + 2) 〕
=〔ª (1) , x (0)(2) , ª (3) , ⋯, ª ( q + 1) , x ( q + 2) 〕
其中 x (0)(2) = x (0)( k 3 )
x (0)( q + 2) = x (0)( k 3 + q)
设环境突变因素的作用是在 k 3 时刻加入到系统中去的 , x
( k 3 + q) 为系统在 k 3 + q 时经综合分析得出的合理期望值或
k 3 + q 时刻的真实值 ,则在 k 3 以后时段系统行为的发展系数
及灰作用量满足
a = - ln ( X (0)( k 3 + q) )X (0)( k 3 )q
u = x (0)( k 3 ) ·ln ( x (0)( k 3 + q)x (0)( k 3 ) )
q[ ( x
(0)
( k 3 + q)
x
(0)
( k 3 ) ) 1/ q - 1 ]
突变预测模型的具体方法请参阅文献[13 ] .
3 结果与分析
3 . 1 气象控制因子选择
由于杉木直径生长受气候条件影响显著 ,根据历
史经验及杉木气候学研究结果[12 ,15 ] ,影响杉木直径生
长的五大主导因子为 : 15cm 和 20cm 年平均地温、全
年日平均气温在 18～27 ℃的日数、5 和 6 月月平均最
小相对湿度及 9 和 10 月月平均降水量. 因此 ,本文以
这 5 个气候因子作为多维时间序列分析模型中的控制
项 ,因为预测某一年直径生长需用上一年以前的资料 ,
因此可以提前一年作出预报 ,直径生长以年轮指数作
为时间序列 ,将年轮指数及气象控制因子列于表 1.
年轮指数及气象资料组建成数据文件 ,连接于用
表 1 直径生长和气象资料
Table 1 Diameter growth and climate data
A B C D E F G
6 23. 3 23. 2 145 37 94. 45 1. 3125
7 23. 3 23. 2 156 43 54. 55 1. 1778
8 22. 6 22. 6 176 26. 5 95. 55 0. 9286
9 23. 3 23. 2 165 31 71. 35 1. 0500
10 22. 7 22. 6 147 34 33. 45 1. 0216
11 23. 3 22. 2 146 30. 5 41. 85 1. 2000
12 23. 3 22. 3 155 42. 5 99. 10 1. 1818
13 21. 8 21. 7 138 36. 5 90. 65 1. 0968
14 22. 5 22. 6 127 29. 5 67. 70 1. 1034
15 22. 5 22. 6 168 36. 0 64. 15 0. 8519
16 22. 1 22. 2 164 48. 0 130. 25 0. 8077
17 22. 2 22. 2 180 38. 5 102. 60 0. 8750
18 22. 0 22. 1 169 28. 0 179. 90 1. 0000
19 21. 8 21. 8 162 32. 5 127. 40 0. 8095
20 22. 6 22. 7 164 30. 0 80. 45 1. 0500
21 22. 4 22. 4 155 30. 0 48. 25 1. 0526
22 22. 2 22. 3 143 22. 5 81. 15 1. 1111
23 22. 7 22. 8 133 27. 5 86. 05 1. 0000
24 22. 5 22. 7 158 29. 5 100. 05 0. 8750
A. 树龄 Age (yr) ,B. 15cm 年均地温 Annual average geotemperature under
15cm( ℃) ,C. 20cm 年均地温 Annual average geotemperature under 20cm
( ℃) ,D. 全年日均气温在 18～20 ℃日数 Days of mean temperature from
18 to 27 in a year (d) ,E. 5 月和 6 月月平均最小相对湿度 Mean minimum
relative humidity in May and J une ( %) , F. 9 月和 10 月月平均降水量
Mean rainfall in September and October ( mm) , G. 年轮指数 Age ring in2
dex.
693 应 用 生 态 学 报 10 卷
BASIC语言编制的 CAR 模型 ,输入相应的参数如资
料年数、F 检验值等 ,即可输出预测模型 ,并作出预
测.
3 . 2 CAR 模型最高阶 n 的 F 检验
由低级开始递增地建立 CAR (1) 、CAR (2) 、CAR
(3) 、CAR (4)模型 ,即建立下列模型 :
Y t =α1 Y t - 1 + ⋯+αn Y t - n + b0 U t + b1 U t - 1 + ⋯
+ bnU t - n +εt (7)
F 检验结果列于表 2 ,CAR (2)与 CAR (3) 模型用 F 检
验法检验得相邻的 F 值等于 1. 1928 ,明显小于 Fα ,因
此 ,CAR (2)模型是合适的.
表 2 F检验结果
Table 2 Results of F2test
参数 Parameter a b F 值 F value
221 0. 2258 0. 7888 605. 29
322 - 0. 0112 0. 9975 1. 1928
3 . 3 CAR 模型真实阶及时滞的 F 检验
虽然 CAR (2) 模型是合适的 ,但其中较多变量参
数接近于零. 现将接近于零的参数剔除并重新建立
CAR (2)模型. 杉木人工林直径生长的年轮指数多维时
间序列模型为 :
Y = 0 . 3915003 Y1 - 0 . 1445521 Y2 +
0. 8720304 X11 - 0 . 9142589 X21 +
0. 07661726 X22 (8)
F1 = 1 . 8605 R = 0 . 9976 F2 = 1 . 1928
只要将相应的气候因子、预测因子实际值代入预
测模型 ,即可对各年份作回归检验或作出下一年的实
际预测预报. 其中 Y1 为上一年的年轮指数 ; Y2 为上
上年的年轮指数 ; X11 为上一年的 15cm 年均地温
( ℃) ; X21为上年的 20cm 年均地温 ( ℃) ; X22为上上年
的 20cm 年均地温 ( ℃) ; F1 为模型参数检验 F 值 ; F2
为相邻阶模型比较 F 值. 从预测模型可以看出 ,杉木
直径生长与 15cm 年均地温以及 20cm 年均地温有关 ,
从而为进一步研究杉木直径生长机制、杉木气候学提
供了统计学生态信息.
3 . 4 突变预测模型的建立
由表 3 可以看出 ,19、20、21 年模拟误差较大 ,这
与全球气候变化频繁有关 ,能否将突变预测模型引入
以降低模拟误差 ,值得探讨. 因此先作地温预测 ,并考
虑突变因素 ,由含突变因素的 GM ( 1 , 1) 模型建立
15cm 年均地温、20cm 年均地温的突变预测模型如下 :
15cm 年均地温 :
X^ (0)( k + 1) = 10136 . 2375 ( e0. 002217298 k -
e0. 002217298 ( k - 1) ) (9)
20cm 年均地温 :
X^ (0)( k + 1) = 10226 . 4875 ( e0. 002207509 k -
e0. 002217298 ( k - 1) ) (10)
其中 ,设环境突变因素的作用是在 15 年加入到系统中
去的.
3 . 5 突变预测模型与多维时间序列模型结合
先利用突变预测模型式 (9) 、(10) 对 15cm 年均地
温、20cm 年均地温进行预测 ,偶后再用多维时间序列
模型建立年轮指数生长的 2 阶 CAR (2)模型 :
Y = 0 . 3040133 Y1 - 0 . 2039253 Y2 +
0. 1078202 X11 - 0 . 2534115 X21 +
0. 1848469 X22 (11)
R = 0 . 9966
3 . 6 回检
利用预测 CAR (2)模型 (8) 和 (11) ,对杉木人工林
直径生长的年轮指数进行回归验证 ,表 3 列出了回检
的结果. 结果表明 ,多维时间序列模型应用于林木生长
表 3 回检结果
Table 3 Simulating results
树龄
Age
(yr)
实际值
Real value
CAR 模拟值 Simulating value of CAR
模型 (8) Model 8
数值 Value △( %)
模型 (11) Model 11
数值 Value △( %)
逐步回归法模拟值
Simulating
value of step
regression method
△( %)
10 1. 0216 1. 1159 9. 23 0. 9405 7. 93 1. 1153 9. 17
11 1. 2000 1. 1585 3. 46 1. 1053 7. 89 1. 1141 7. 16
12 1. 1818 1. 2034 1. 83 1. 2205 3. 27 0. 9839 16. 75
13 1. 0968 1. 0484 4. 41 1. 0793 1. 59 0. 9821 10. 46
14 1. 1034 1. 1380 3. 14 1. 0659 3. 39 0. 9224 16. 40
15 0. 8519 0. 8945 5. 00 0. 8218 3. 53 1. 0553 23. 88
16 0. 8077 0. 8640 6. 97 0. 9030 11. 79 0. 8947 10. 77
17 0. 8750 0. 8999 2. 85 0. 9429 7. 76 1. 1440 30. 74
18 1. 0000 0. 9892 1. 08 0. 9743 2. 57 1. 1469 14. 69
19 0. 8095 0. 9455 16. 80 1. 0006 23. 61 1. 0916 34. 85
20 1. 0500 0. 9450 10. 00 0. 9191 12. 46 1. 1163 6. 31
21 1. 0526 0. 9185 12. 74 1. 0331 1. 85 1. 1080 5. 26
22 1. 1111 1. 0536 5. 18 0. 9868 11. 18 1. 0554 5. 02
23 1. 0000 0. 9702 2. 98 1. 0059 0. 59 1. 0002 0. 02
24 0. 8750 0. 8894 1. 65 0. 9622 9. 96 1. 0219 16. 78
精度 Precision ─ ─ 94. 18 ─ 92. 75 ─ 86. 12
△= │模拟值 - 真实值│/ 真实值│Simulating value - Real value │/ Real value.
7934 期 吴承祯等 :林木生长的多维时间序列分析
序列分析是可行的 ,其拟合精度较高. 与用逐步回归法
相比较[12 ] ,其回归优度可依次提高 4. 8 和 4. 2 个百分
点 ,精度提高 8. 06 和 6. 63 个百分点. 就模型 (8) 和模
型 (11)而言 ,引进自变量的突变预测模型后 ,CAR (2)
模型回检精度比原模型差 ,精度降低了 1. 43 个百分点
且误差变幅更大 ,可以认为以式 (8)作为杉木人工林直
径生长的年轮指数多维时间序列模式模型更合理、可
靠 ,因此杉木人工林直径生长的年轮指数多维时间序
列模型为 :
Y = 0 . 3915003 Y1 - 0 . 1445521 Y2 +
0. 8720304 X11 - 0 . 9142589 X21 +
0. 07661726 X22 (12)
4 讨 论
根据杉木直径生长年轮指数的 24 年资料以及相
应各年的气候组合因子资料 ,应用逐步回归分析方法
筛选对杉木直径生长起主导作用的气候因子为 :15cm
年均地温、20cm 年均地温、全年日平均气温在 18～
27 ℃的日数、5 和 6 月月平均最小相对湿度及 9 月和
10 月月平均降水量[12 ,15 ] . 本文以这 5 个主导气候因
子为多维时间序列模型的控制因子 ,建立了杉木直径
生长的 CAR 模型 ,从而提前一年对杉木直径生长作出
科学预测. 但由于全球气候异常变化 ,气候因子出现相
应异常值 ,从而使某些年预报误差较大 ,因此有必要先
对控制因子进行预测 ,并考虑突变因素 ,再进行建模.
但引入突变预测模型后效果并不理想 ,有待深入研究.
杉木直径生长的多维时间序列 CAR 模型的建模
关键在于选择控制因子 ,由于杉木气候学对杉木生长
系统的内部机制做了深入的研究 ,寻找到了合适的因
子 ,因此所建立的多维时间序列 CAR 模型在低阶时模
型效果就很好 ,故可以比较有把握地探知林木直径生
长的未来 ,从而为林木生长预测提供新方法 ,也为树木
气候学研究引入新内容、提供新信息.
CAR 模型结合了一维时间序列分析和回归分析
两数数理统计方法的优点 ,不仅考虑到事物发展的自
身运动规律 ,也顾及了环境因子的作用 ,因而可以认为
是一个值得探讨和应用的方法.
致谢 本文在修改过程中得到北京林业大学袁嘉祖教授的指
导和帮助 ,在此深表感谢 !
参考文献
1 王振中、林孔勋、范怀 . 1990. 小白菜花叶病时间序列分析预测
法. 植物保护学报 ,19 (3) :269～273.
2 云南林学院主编. 1979. 气象学. 北京 :农业出版社.
3 邓自立. 1992. 动态系统分析及其应用. 沈阳 :辽宁科学技术出版
社. 76～191.
4 兰 斌. 1990. 地位指数时间序列建模及其应用研究. 福建林学院
学报 ,10 (2) :146～151.
5 李鸿杰. 1993. 土壤水、盐、入渗变异特性及其相互关系的空间序列
分析. 土壤学报 ,30 (1) :60～68.
6 周立阳. 1995. 多维时间序列分析在稻纵卷叶螟长期预测预报上的
应用. 植物保护学报 ,22 (1) :1～6.
7 张荷观. 1986. 时间序列的 ARIMA 模型在预测树木生长量中的应
用. 林业科学 ,22 (1) :94～100.
8 邱学清、江希钿. 1989. AR( P) 模型在预测黑荆树林分月总生长量
中的应用. 福建林学院学报 ,9 (1) :43～49.
9 余新晓. 1990. 降雨侵蚀力指数的时间序列分析. 北京林业大学学
报 ,12 (2) :55～62.
10 张素芬、屠乃斌. 1992. 油松毛虫种群动态的 ARMA ( P ,q) 模型. 北
京林业大学学报 ,14 (1) :93～97.
11 张明芝、汪寅虎、柯福源. 1990. 有机肥施用量预测研究. 土壤通报 ,
21 (4) :185～188.
12 吴 敬、洪 伟. 1983. 杉木气候学研究. 气象科技 , (3) :47～54.
13 袁嘉祖. 1992. 灰色系统及其应用. 北京 :科学出版社. 81～83.
14 姚晓红. 1990. 林分生长数据的时间分析探讨. 北京林业大学学报 ,
12 (4) :10～16.
15 洪 伟、林思祖编著. 1993. 计量林学研究. 成都 :电子科技大学出
版社. 175～178.
16 程极益. 1992. 作物病虫害数理统计预报. 北京 :农业出版社. 75～
182.
17 G. W. Cox. 1979. 普通生态学实验手册. 北京 :科学出版社.
18 Box , G. E. P. 1970. Time series analysis : Forcasting and control.
Holden2Day , San Francisco ,USA.
作者简介 吴承祯 ,男 ,29 岁 ,副教授 ,在职博士 ,主要从事数量
生态和林业系统工程、环境科学等领域研究 ,发表学术论文 70
余篇.
893 应 用 生 态 学 报 10 卷