CN108776850A - 一种森林蓄积量的精确预估方法 - Google Patents

Abstract

本发明公开了一种森林蓄积量的精确预估方法，按如下步骤进行：一是收集对象地域的时间序列成对数据；二是确定自回归滑动平均模型的阶；三是确定向量矩阵模型；四是求参数估计量；五是新息估计；六是扩展参数估计；七是新息再估计；八是扩展参数再估计；九是扩展参数最终估计；十是蓄积量的估计与预测。本发明具有五方面的优点：能综合实测不能得到的影响蓄积量因素；既能描述蓄积量的自相关性，又能反映蓄积量与有林地面积之间的回归关系；能反映树木生长的非线性；有完备的理论基础，能克服经验模型的缺陷；能大大消弱新息序列对模型参数估计的影响，提高模型估计与预测精度。

Description

一种森林蓄积量的精确预估方法

技术领域

本发明涉及一种森林蓄积量估计与预测方法，特别是一种基于自回归滑动平均模型与扩展最小二乘的森林蓄积量估计与预测方法。

背景技术

森林蓄积量是区域自然资源的重要组成部分,也是森林资源调查和监测及预测的最重要的指标之一。高效、快速、准确地估计与预测森林蓄积量，是评估森林碳收支的基础,还能为森林资源可持续发展和区域林业发展战略制定提供决策依据。快速、精确地估计与预测森林蓄积量，属于面向决策的科学范畴，采用何种估计与预测方法很重要。

目前有关森林蓄积的估计与预测方法主要有四种：一是神经网络，该方法的优点是适于模拟复杂系统,且具有自学习、联想存贮和高速寻找优化解的功能，主要缺陷是收敛速度慢且难以确定隐层及隐节点的个数，估算精度受其影响；二是灰色系统理论，该方法首先对原始非平稳时间序列难以取得满意的预测效果，其次已知条件缺乏严格的理论依据，故估算结果不一定是最佳预测；三是k近邻法(K-NN)，k-NN方法曾用于对意大利两个地区森林蓄积量进行研究，并且研究表明估计与预测结果受遥感数据光谱类型和输入辅助变量、所测多维距离的类型、近邻样本个数等因素的影响；四是3S技术的应用，这种方法的优点是为大尺度森林蓄积量估计与预测提供了一条快捷、经济、方便和可靠的途径，缺陷是遥感数据有其局限性，特别是在多云雨和雾的热带和亚热带地区,可见光和红外遥感受到了很大限制，并且蓄积量估计与预测变量一般选取如林龄、郁闭度、土壤、土壤厚度、优势树种等因素，但这些因素很难通过遥感、GIS获取，此外蓄积量估计与预测模型大多选取多项式模型，由于森林分布、类型及生长状况等相差较大，再加之多项式方程是经验模型，带有很大的人为主观性，因此多项式模型不一定适合于森林蓄积量的估计与预测。

现有的森林蓄积量预测和预估方法在计算过程中受噪声污染的影响较大，结果往往有不小的偏差。因此，设计一种减弱甚至消除噪声污染，结果精确的，可以充分利用森林资源调查资料的信息，并且客观反映森林资源总规模及其丰富程度，进而能衡量森林生态环境的优劣及给区域林业发展规划提供可靠的数据支撑的森林蓄积量估计与预测方法是本行业技术人员继续解决的问题。

发明内容

本发明的目的在于，提供森林蓄积量的无偏预估方法。本发明有效消弱了有色噪声对参数估计的影响，提高森林蓄积量的估计与预测精度。同时可以充分利用森林资源调查资料的信息，能客观并且相对精确地反映森林生长情况，进而可以衡量森林生态环境的优劣及给区域林业发展规划提供可靠的数据支撑。

本发明的技术方案：一种森林蓄积量的精确预估方法，包括以下步骤：

1)、收集对象地域最近连续的p期有林地面积与森林蓄积量实测时间序列成对数据，且p＞3n+1，n为公式(1)中模型的阶；

2)、根据步骤1)的数据，基于自回归滑动平均模型与扩展最小二乘算法，计算森林蓄积量，具体步骤如下：

a)森林蓄积量时间序列模型：

设森林蓄积量时间序列模型为：

其中：n为模型的阶；k、j均为常数且k>j；

y(k)为系统的输出，代表第k期的森林蓄积量；

y(k-j)为系统的历史输出，代表第k-j期的森林蓄积量；

u(k-j)为系统可确定性输入，代表第k-j期的有林地面积；

ε(k),ε(k-1),…,ε(k-n)为新息序列，即该序列的均值为零且互不相关；

a_j为自回归系数，b_j为输入传递系数，c_j为新息系数,即a_j、b_j与c_j为森林蓄积量时间序列模型的参数；

b)确定模型阶n：

令θ＝(a₁,a₂,...,a_n,b₀,b₁...,b_n)，代表公式(1)式中森林蓄积量与有林地面积的参数向量；再令c₁＝c₂＝…＝c_n＝0，用不同阶次的模型进行最小二乘拟合，计算拟合优度函数

其中是在某一阶次下模型参数θ的估计量；

当模型阶次增加时，的值减小，在显著减小终止时，对应模型的阶次即为合适的模型阶次；

c)写出向量矩阵模型：

令β＝(a₁,a₂,...,a_n,b₀,b₁...,b_n,c₁,c₂,…,c_n)；其中：为森林蓄积量的n期历史输出信息；

为有林地的n+1期输入信息；

为n期新息序列；

β为公式(1)所示模型的参数向量；

则公式(1)的另种表示形式为：

y(k)＝x^T(k)β^T+ε(k)； (2)

当k＝n+1,n+2,...,n+p时，(2)式的向量矩阵模型为：

Y＝Xβ^T+ε； (3)

其中：Y^T＝(y(n+1),y(n+2),...,y(n+p))，

ε^T＝(ε(n+1),ε(n+2),...,ε(n+p))，

d)求参数θ的估计量即扩展最小二乘首先需求的量：

令c₁＝c₂＝…＝c_n＝0，则公式(1)变为：

再令则公式(4)的矩阵形式为：

当k＝n+1,n+2,...,n+p时，

其中：

e)求新息的估计值：

由公式(5)知，则新息的估计值为：

f)求参数β估计值

由公式(7)得

其中为第k-1期、第k-2期、…、第k-n期新息的估计值；

则

由公式(3)得为：

其中：为的转置；

g)再求新息的估计值：

由公式(2)得ε(k)＝y(k)-x^T(k)β^T，根据步骤f)得则新息序列的估计值为：

h)再求参数β的估计值

由步骤g)得进而得及此时为：

i)参数β最终估计值：

用步骤h)中β的新估计值，代替步骤g)中β的估计值；再次进行步骤g)和步骤h)，重新计算及后再得β的最新估计值；

如此循环，直到β的估计值稳定为止；

j)森林蓄积量的估计与预测：

将步骤i)最终得到的代入公式(1)，得y(k)估计值或预测值，即得对象地域第k期森林蓄积量的估计值或预测值。

与现有技术相比，本申请从如下两方面着手来精确估计与预测森林蓄积量：

一是时间序列分析。这是一种广泛应用的数据分析方法,它主要用于描述和探索现象随时间发展变化的数量规律。时间序列分析的一项重要内容就是根据过去已有的数据来估计现状与预测未来。由于影响森林蓄积的因素错综复杂且有些影响因素的数据资料无法得到，这时以时间t综合替代这些因素的时间序列分析方法就能达到估计与预测的目的。自回归滑动平均模型即ARMA模型是依据森林生长自身的量测因子时间序列建立起来的，而不用管生态系统内部的运行情况，因此是一种解决非确定性系统的较好模型,它不仅能对系统作出状态估计,也能对系统内部的运行状况进行趋势预测。再加之，森林蓄积量是与时间及有林地面积有关的变量，因此估计与预测森林蓄积量的首选模型是时间序列分析的ARMA模型。因每期资源调查的有林地面积与森林蓄积量是成对的采样数据，故选用ARMA(n,n)模型；

二是用扩展最小二乘辨识ARMA模型各参数。ARMA模型最常用的参数辨识方法是普通最小二乘算法，普通最小二乘算法只有在ARMA模型残差项为白噪声时对蓄积量的估计才是无偏的。但由于受很多自然与非自然因素的随机干扰，使森林生态系统内部特征具有随机性，这样ARMA模型的输入端必然会受到噪声的污染，其残差项在一般情况下为有色噪声。故用普通最小二乘算法估计与预测森林蓄积量的结果是有偏的，而扩展最小二乘把有色噪声看作新息序列，并估算出新息序列的参数，所以得到的结果比普通最小二乘精确，即通过扩大观测矩阵，消弱新息序列对模型参数估计的影响，从而提高蓄积量的估计与预测精度。

由于森林资源是不断变化的，考虑了随机干扰给预测模型本身带来的影响。因此本申请用随机变量来描述森林变化，使结果更加精确。自回归滑动平均模型也可以看作是系统的差分方程模型，而差分方程可以描述一些符合树木生长的非线性模型。同时，本申请模型具有完备的理论基础，可以克服有些估算方法选用经验模型的缺陷。

综上，发明有效消弱了有色噪声对参数估计的影响，提高森林蓄积量的估计与预测精度。同时可以充分利用森林资源调查资料的信息，能客观并且相对精确地反映森林生长情况，进而可以衡量森林生态环境的优劣及给区域林业发展规划提供可靠的数据支撑。

具体实施方式

下面结合实施例对本发明作进一步说明，但并不作为对本发明限制的依据。

实施例：一种森林蓄积量的精确预估方法，包括以下步骤：

1)、收集对象地域最近连续的p期有林地面积与森林蓄积量实测时间序列成对数据，且p＞3n+1，n为公式(1)中模型的阶。在收集数据时，数据越多越好，也即p值越大越好。而n的值又是可预估的，经过拟合通常取3或4。因此在收集数据时，基本都能满足p＞3n+1要求。

2)、根据步骤1)的数据，基于自回归滑动平均模型与扩展最小二乘算法，计算森林蓄积量，具体步骤如下：

a)森林蓄积量时间序列模型：

设森林蓄积量时间序列模型为：

其中：n为模型的阶；k、j均为常数且k>j；

y(k)为系统的输出，代表第k期的森林蓄积量；

y(k-j)为系统的历史输出，代表第k-j期的森林蓄积量；

u(k-j)为系统可确定性输入，代表第k-j期的有林地面积；

ε(k),ε(k-1),…,ε(k-n)为新息序列，即该序列的均值为零且互不相关；

a_j为自回归系数，b_j为输入传递系数，c_j为新息系数,即a_j、b_j与c_j为森林蓄积量时间序列模型的参数；

b)确定模型阶n：

令θ＝(a₁,a₂,...,a_n,b₀,b₁...,b_n)，代表公式(1)式中森林蓄积量与有林地面积的参数向量；再令c₁＝c₂＝...＝c_n＝0，用不同阶次的模型进行最小二乘拟合，计算拟合优度函数

其中是在某一阶次下模型参数θ的估计量；

当模型阶次增加时，的值减小，在显著减小终止时，对应模型的阶次即为合适的模型阶次；

c)写出向量矩阵模型：

令β＝(a₁,a₂,...,a_n,b₀,b₁...,b_n,c₁,c₂,…,c_n)；其中：为森林蓄积量的n期历史输出信息；

为有林地的n+1期输入信息；

为n期新息序列；

β为公式(1)所示模型的参数向量；

则公式(1)的另种表示形式为：

y(k)＝x^T(k)β^T+ε(k)； (2)

当k＝n+1,n+2,...,n+p时，(2)式的向量矩阵模型为：

Y＝Xβ^T+ε； (3)

其中：Y^T＝(y(n+1),y(n+2),...,y(n+p))，

ε^T＝(ε(n+1),ε(n+2),...,ε(n+p))，

d)求参数θ的估计量即扩展最小二乘首先需求的量：

令c₁＝c₂＝…＝c_n＝0，则公式(1)变为：

再令则公式(4)的矩阵形式为：

当k＝n+1,n+2,...,n+p时，

其中：

e)求新息的估计值：

由公式(5)知，则新息的估计值为：

f)求参数β估计值

由公式(7)得

其中为第k-1期、第k-2期、…、第k-n期新息的估计值；

则

由公式(3)得为：

其中：为的转置；

g)再求新息的估计值：

由公式(2)得ε(k)＝y(k)-x^T(k)β^T，根据步骤f)得则新息序列的估计值为：

h)再求参数β的估计值

由步骤g)得进而得及此时为：

i)参数β最终估计值：

用步骤h)中β的新估计值，代替步骤g)中β的估计值；再次进行步骤g)和步骤h)，重新计算及后再得β的最新估计值；

如此循环，直到β的估计值稳定为止；

j)森林蓄积量的估计与预测：

将步骤i)最终得到的代入公式(1)，得y(k)估计值或预测值，即得对象地域第k期森林蓄积量的估计值或预测值。

应用实例如下：

用1986-1999年福州市永泰县城峰镇有林地面积、森林蓄积量数据,进行分析说明，其具体数值见表1：

表1城峰镇1986-1999年有林地面积(1000hm²)、森林蓄积量(10000m³)

一、不同方法的计算结果

经确定ARMA模型阶次为3。

假设时间退回到1997年，则1986-1997年的有林地面积与森林蓄积量是收集的实测时间序列成对数据，而1998、1999年有林地面积0.523、0.534万hm²可由该镇林业规划事先确定的，这样对1989-1997年城峰镇森林蓄积量，用4种方法得到的是估计值，而对1998、1999年的蓄积量，4种方法得到的则是预测值。

普通最小二乘参数为：

(a₁,a₂,a₃,b₀,b₁,b₂,b₃)＝(-0.1012,-0.3477,-0.3915,0.4300,0.3026,1.1027,-1.3297)

则扩展最小二乘经70次迭代后的参数为：

(a₁,a₂,a₃,b₀,b₁,b₂,b₃,c₁,c₂,c₃)

＝(-0.2051,-0.2208,0.3043,-0.6731,-1.5777,14.2755,-8.8467,-0.0423,-0.08,0.0002)

采用一个隐含层的前馈BP神经网络对蓄积量的估计与预测结果见表2，其中隐含层有5个神经元，输出层有1个神经元；

由于森林蓄积量是与有林地面积有关的变量，故采用GM(1,2)的灰色模型进行估计与预测，经计算得GM(1,2)模型为

解此微分方程得如下函数：

城峰镇森林蓄积量的实际值及每种方法的估计值与预测值见表2。

表2四种方法的蓄积量估计与预测结果及与实测值的对比

二、结果分析

普通最小二乘算法不考虑新息序列，所以也不能得出新息序列的参数，因此得到的估计值、预测值是有偏的，而对1998与1999年森林蓄积量的预测偏差较大。

扩展最小二乘法考虑了新息序列对蓄积量估计与预测的影响，因此大大消弱了新息序列对参数的估计的影响，从而提高了估计与预测精度。

GM(1,2)模型的估计值、预测值的残差方差及对蓄积量的预测可以看出，GM(1,2)模型的估计与预测效果最差，主要是由于GM(1,2)是用的指数曲线去拟合因此其对蓄积量的估计与预测完全依赖于有林地面积，其实蓄积量是与时间有关的量即每期蓄积量之间也存在着一定的自回归关系，由于灰色模型没考虑蓄积量之间的自相关性，再加之，在残差序列是有色噪声时采用普通最小二乘算法进行参数辨识，因此其估计与预测精度最低。

BP神经网络的估计值、预测值的残差方差，及对1989-1997、1998、1999年森林蓄积量的估计与预测可以看出，该方法对信息的识别能力比较差，主要是由于：1)BP神经网络对初始权重非常敏感，极易收敛于局部极小，2)网络隐含节点数的确定至今只有一些经验公式，而没有任何理论上的指导。

Claims (1)

1.一种森林蓄积量的精确预估方法，其特征在于：包括以下步骤：

1)、收集对象地域最近连续的p期有林地面积与森林蓄积量实测时间序列成对数据，且p＞3n+1，n为公式(1)中模型的阶；

2)、根据步骤1)的数据，基于自回归滑动平均模型与扩展最小二乘算法，计算森林蓄积量，具体步骤如下：

a)森林蓄积量时间序列模型：

设森林蓄积量时间序列模型为：

其中：n为模型的阶；k、j均为常数且k>j；

y(k)为系统的输出，代表第k期的森林蓄积量；

y(k-j)为系统的历史输出，代表第k-j期的森林蓄积量；

u(k-j)为系统可确定性输入，代表第k-j期的有林地面积；

ε(k),ε(k-1),…,ε(k-n)为新息序列，即该序列的均值为零且互不相关；

a_j为自回归系数，b_j为输入传递系数，c_j为新息系数,即a_j、b_j与c_j为森林蓄积量时间序列模型的参数；

b)确定模型阶n：

令θ＝(a₁,a₂,...,a_n,b₀,b₁...,b_n)，代表公式(1)式中森林蓄积量与有林地面积的参数向量；再令c₁＝c₂＝…＝c_n＝0，用不同阶次的模型进行最小二乘拟合，计算拟合优度函数

其中是在某一阶次下模型参数θ的估计量；

当模型阶次增加时，的值减小，在显著减小终止时，对应模型的阶次即为合适的模型阶次；

c)写出向量矩阵模型：

令β＝(a₁,a₂,...,a_n,b₀,b₁...,b_n,c₁,c₂,…,c_n)；其中：为森林蓄积量的n期历史输出信息；

为有林地的n+1期输入信息；

为n期新息序列；

β为公式(1)所示模型的参数向量；

则公式(1)的另种表示形式为：

y(k)＝x^T(k)β^T+ε(k)； (2)

当k＝n+1,n+2,...,n+p时，(2)式的向量矩阵模型为：

Y＝Xβ^T+ε； (3)

其中：Y^T＝(y(n+1),y(n+2),...,y(n+p))，

ε^T＝(ε(n+1),ε(n+2),...,ε(n+p))，

d)求参数θ的估计量即扩展最小二乘首先需求的量：

令c₁＝c₂＝…＝c_n＝0，则公式(1)变为：

再令则公式(4)的矩阵形式为：

当k＝n+1,n+2,...,n+p时，

其中：

e)求新息的估计值：

由公式(5)知，则新息的估计值为：

f)求参数β估计值

由公式(7)得

其中为第k-1期、第k-2期、…、第k-n期新息的估计值；

则

由公式(3)得为：

其中：为的转置；

g)再求新息的估计值：

由公式(2)得ε(k)＝y(k)-x^T(k)β^T，根据步骤f)得则新息序列的估计值为：

h)再求参数β的估计值

由步骤g)得进而得及此时为：

i)参数β最终估计值：

用步骤h)中β的新估计值，代替步骤g)中β的估计值；再次进行步骤g)和步骤h)，重新计算及后再得β的最新估计值；

如此循环，直到β的估计值稳定为止；

j)森林蓄积量的估计与预测：

将步骤i)最终得到的代入公式(1)，得y(k)估计值或预测值，即得对象地域第k期森林蓄积量的估计值或预测值。