森林蓄积量的无偏预估方法与流程

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

背景技术

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

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

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

技术实现要素：

本发明的目的在于，提供森林蓄积量的无偏预估方法。本发明有效减弱甚至消除了噪声污染的影响，结果更加精确，同时可以充分利用森林资源调查资料的信息，并且能客观反映森林资源总规模及其丰富程度，进而可以衡量森林生态环境的优劣及给区域林业发展规划提供可靠的数据支撑。

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

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

2)通过步骤1)收集的数据，基于自回归滑动平均模型与广义最小二乘算法，预估对象地域的森林蓄积量，具体如下：

a、森林蓄积量时间序列模型与噪声滞后算子多项式：

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

其中：n为模型的阶；k、j均为常数且k＞j，y(k)为系统的输出，代表第k期的森林蓄积量；y(k-j)为系统的历史输出，代表第k-j期的森林蓄积量；u(k-j)为系统可确定性输入，代表第k-j期的有林地面积；aj为自回归系数，bj为输入传递系数，即aj与bj为森林蓄积量时间序列模型的参数；ξ(k)为第k期的残差，即模型的有色噪声；

设噪声滞后算子多项式为：

其中：z^-i(i＝1,2,...,n)为滞后算子，a1,a2,…,an,b0,b1,b2,…,bn为参数；

则公式(1)可表示为：

a(z^-1)y(k)＝b(z^-1)u(k)+ξ(k)(3)

b、确定模型阶n：

令θ＝(a1,a2,…,an,b0,b1...,bn)，代表模型参数向量；

用不同阶次的模型进行最小二乘拟合，计算拟合优度

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

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

c、计算参数的一次估计量，即广义最小二乘算法需先求的量：

设白化滤波器即白化函数为：

f(z^-1)＝1+f1z^-1+...+fmz^-m；(4)

其中：z^-j(j＝1,2,...,m)为滞后算子，f1,f2,…,fm为参数；

有色噪声ξ(k)与白噪声w(k)的转化关系如下：

其中：f(z^-1)为白化滤波器；

令公式(5)中的f(z^-1)＝1，即噪声模型参数向量f为m维零向量，其中f＝(f1,f2,...,fm)，则模型a(z^-1)y(k)＝b(z^-1)u(k)+ξ(k)或的残差便退化为白噪声，即ξ(k)＝w(k)；

令

其中：为n期历史输出信息，为n+1期输入信息；

当得到p对采样值{x^t(k),y(k)},k＝n+1,n+2,...,n+p时，根据y(k)＝x^t(k)θ^t+ξ(k)，得出向量矩阵模型：

y＝xθ^t+ξ；(6)

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

按普通最小二乘法求得θ^t的一次近似估计量：

d、求噪声模型参数向量：

令φ^t(k)代表森林生态系统随机干扰下的有色噪声序列；利用及各期的成对样值，由y(k)＝x^t(k)θ^t+ξ(k)算出模型各期有色噪声ξ(k)的估计值：

将得到的分别代入得到：

根据向量矩阵模型：

其中：

利用利用普通最小二乘法求得f^t的估计量：

其中：为的转置；

e、计算滤波信号：

将公式(5)代入a(z^-1)y(k)＝b(z^-1)u(k)+ξ(k)，得

a(z^-1)f(z^-1)y(k)＝b(z^-1)f(z^-1)u(k)+w(k)；

利用步骤d得到的计算输入滤波信号：

与输出滤波信号：

以及输入输出滤波信号：

f、求参数的二次估计量：

根据步骤d求得再根据步骤e，求得与

用与分别代替公式f(z^-1)y(k)＝f(z^-1)x^t(k)θ^t+w(k)中的y(k)与x^t(k)，得到输出滤波信号与输入输出滤波信号间的普通线性回归模型：

当k＝n+1,n+2,…,n+p时，得向量矩阵模型：

其中：

利用普通最小二乘法求得θ^t的二次近似估计量：

g、计算参数最终估计量：

用取代根据步骤d，计算

根据步骤f，计算得到参数的三次近似估计量

如此反复迭代，直到其中：i为的第i次迭代，n为的迭代总次数；

此时，θ^t的估计量再将代入

得y(k)的估计值或预测值，即得对象地域第k期森林蓄积量的估计值或预测值。

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

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

二是用广义最小二乘算法辨识arma模型各参数。arma模型最常用的参数辨识方法是普通最小二乘算法，然而普通最小二乘算法只有在arma模型残差项为白噪声时对蓄积量的估计才是无偏的。由于受很多自然与非自然因素的随机干扰，森林生态系统内部特征具有随机性，这样arma模型的输入端必然会受到噪声的污染，其残差项在一般情况下为有色噪声。故用普通最小二乘算法估计与预测森林蓄积量的结果是有偏的。而广义最小二乘算法把有色噪声滤波成白噪声，即把ξ(k)转化为w(k)，进而得到森林蓄积量的无偏估计值，提高森林蓄积量的估计与预测精度。

本申请的模型能充分反映森林自身的生长规律并且相对精确。森林生态系统是随机系统，其蓄积量与很多因素有关，且下期的蓄积是前几期蓄积的积累或消耗，因此蓄积之间存在一定的自相关关系，故用时间t综合了实测不能得到的影响蓄积因素，从而建立的与有林地面积有关的arma蓄积模型能反映蓄积的变化规律。常用森林蓄积预估方法一般都是采用确定性定常系统的方法预测森林动态，这与森林资源变化的本身特征不甚相符，而本申请考虑了随机干扰给预测模型本身带来的影响，用随机变量来描述森林变化，使结果更加精确。

森林的生长序列是一个非平稳序列,因此要精确描述这一过程是比较复杂和有实际困难的。而自回归滑动平均模型是以森林自身状态特征为主导信息的动态模型,它能精确描述非平稳状态下森林生长的的积累或消耗过程。同时，该模型具有完备的理论基础，可以克服有些估算方法选用经验模型的缺陷。自回归滑动平均模型也可以看作是系统的差分方程模型，而差分方程可以描述一些符合树木生长的非线性模型，如令y(t)＝e^t，则再如令y(t)＝t^a(a是整数)，则其中为参数，从而(t+1)^a也可以用差分方程来表示。模型结合的广义最小二乘能得到参数的无偏估计，可以克服普通最小二乘估计有偏的缺陷。

综上，本发明有效减弱甚至消除了噪声污染的影响，结果更加精确，同时可以充分利用森林资源调查资料的信息，并且能客观反映森林资源总规模及其丰富程度，进而可以衡量森林生态环境的优劣及给区域林业发展规划提供可靠的数据支撑。

具体实施方式

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

实施例：森林蓄积量的无偏预估方法，包括以下步骤：

1)收集对象地域最近连续的p期有林地面积与森林蓄积量实测时间序列成对数据，且p＞2n+1，n为公式(1)中模型的阶。在收集成对数据的时候，期数p越多越好，加之n的值在可预估的范围(如n＝2或n＝3等)，因此收集的期数p基本都能满足p＞2n+1。

2)通过步骤1)收集的数据，基于自回归滑动平均模型与广义最小二乘算法，预估对象地域的森林蓄积量，具体如下：

a、森林蓄积量时间序列模型与噪声滞后算子多项式：

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

其中：n为模型的阶；k、j均为常数且k＞j，y(k)为系统的输出，代表第k期的森林蓄积量；y(k-j)为系统的历史输出，代表第k-j期的森林蓄积量；u(k-j)为系统可确定性输入，代表第k-j期的有林地面积；aj为自回归系数，bj为输入传递系数，即aj与bj为森林蓄积量时间序列模型的参数；ξ(k)为第k期的残差，即模型的有色噪声；

设噪声滞后算子多项式为：

其中：z^-i(i＝1,2,...,n)为滞后算子，a1,a2,…,an,b0,b1,b2,…,bn为参数；

则公式(1)可表示为：

a(z^-1)y(k)＝b(z^-1)u(k)+ξ(k)(3)

b、确定模型阶n：

令θ＝(a1,a2,...,an,b0,b1...,bn)，代表模型参数向量；

用不同阶次的模型进行最小二乘拟合，计算拟合优度

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

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

c、计算参数的一次估计量，即广义最小二乘算法需先求的量：

设白化滤波器即白化函数为：

f(z^-1)＝1+f1z^-1+…+fmz^-m；(4)

其中：z^-j(j＝1,2,...,m)为滞后算子，f1,f2,…,fm为参数；

有色噪声ξ(k)与白噪声w(k)的转化关系如下：

其中：f(z^-1)为白化滤波器；

令公式(5)中的f(z^-1)＝1，即噪声模型参数向量f为m维零向量，其中f＝(f1,f2,...,fm)，则模型a(z^-1)y(k)＝b(z^-1)u(k)+ξ(k)或的残差便退化为白噪声，即ξ(k)＝w(k)；

令

其中：为n期历史输出信息，为n+1期输入信息；

当得到p对采样值时，根据y(k)＝x^t(k)θ^t+ξ(k)，得出向量矩阵模型：

y＝xθ^t+ξ；(6)

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

按普通最小二乘法求得θ^t的一次近似估计量：

d、求噪声模型参数向量：

令φ^t(k)代表森林生态系统随机干扰下的有色噪声序列；利用及各期的成对样值，由y(k)＝x^t(k)θ^t+ξ(k)算出模型各期有色噪声ξ(k)的估计值：

将得到的分别代入得到：

根据向量矩阵模型：

其中：

利用利用普通最小二乘法求得f^t的估计量：

其中：为的转置；

e、计算滤波信号：

将公式(5)代入a(z^-1)y(k)＝b(z^-1)u(k)+ξ(k)，得

a(z^-1)f(z^-1)y(k)＝b(z^-1)f(z^-1)u(k)+w(k)；

利用步骤d得到的计算输入滤波信号：

与输出滤波信号：

以及输入输出滤波信号：

f、求参数的二次估计量：

根据步骤d求得再根据步骤e，求得与

用与分别代替公式f(z^-1)y(k)＝f(z^-1)x^t(k)θ^t+w(k)中的y(k)与x^t(k)，得到输出滤波信号与输入输出滤波信号间的普通线性回归模型：

当k＝n+1,n+2,…,n+p时，得向量矩阵模型：

其中：

利用普通最小二乘法求得θ^t的二次近似估计量：

g、计算参数最终估计量：

用取代根据步骤d，计算

根据步骤f，计算得到参数的三次近似估计量

如此反复迭代，直到其中：i为的第i次迭代，n为的迭代总次数；

此时，θ^t的估计量再将代入

得y(k)的估计值或预测值，即得对象地域第k期森林蓄积量的估计值或预测值。

从以上步骤可以看出，广义最小二乘算法的如下特点：1)在估计arma模型参数向量θ^t的同时，还须估计噪声模型参数向量f^t；2)是一种逐次逼近算法，具体计算步骤首先假设f^t已知，用普通最小二乘法求出θ^t的估计量然后用给出f^t的估计量再用去修正如此反复迭代，直到为止；3)在每一步迭代过程中，都采用普通最小二乘法来实现。

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

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

一、不同方法的计算结果

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

令m＝2(m见步骤c公式(4))，则普通最小二乘参数为：

(a1,a2,a3,a4,b0,b1,b2,b3,b4)＝(-0.3332,-0.1201,-0.5102,0.6819,1.2428,-0.0022,4.1678,0.2928,-3.1334)(f1,f2)＝(0.0006,-0.0024)

广义最小二乘参数为：

(a1,a2,a3,a4,b0,b1,b2,b3,b4)＝(0.0816,-0.1260,-1.1647,0.1601,-2.8195,-1.3008,4.7038,2.8720,-3.7722)(f1,f2)＝(0,0)

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

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

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

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

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

二、结果分析

从普通最小二乘算法得出的噪声参数向量f与估计值、预测值的残差方差可以得出：1)残差序列ξ(k)是有色噪声，2)普通最小二乘对参数的估计是有偏的，因此arma模型与普通最小二乘对蓄积量的估计结果是有偏的，而且对1998与1999年森林蓄积量的预测偏差较大。

广义最小二乘算法得出的噪声参数向量f与估计值、预测值的残差方差可以得出：1)有色噪声ξ(k)已经转换为白噪声w(k)，2)广义最小二乘对参数的估计是无偏的，因此arma模型与广义最小二乘对蓄积量的估计结果是无偏的，而对1998与1999年森林蓄积量的预测精度最高。

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

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