CN104615875B - 一种稳健遥感单木冠幅与林木直径回归方法 - Google Patents

Abstract

一种稳健遥感单木冠幅与林木直径回归方法，属于计算机程序及测量技术领域。含有以下步骤：冠幅与林木直径相关数据获取；稳健回归模型的构建；单木冠幅与林木直径的稳健回归模型参数的迭代确定。本发明在典型的东北中温带针阔叶混交林区，在高空间分辨率的单木树冠影像的提取基础上，为克服传统最小二乘法对“异常值点”的数据过于敏感的“缺陷”，采用Huber提出M‑估计法，成功地构造了遥感影像单木冠幅与林木直径稳健回归模型。根据M‑估计法引进目标函数的加权因子分布，从而改变传统的最小二乘法各样本(包括异常点)权重等同的“缺陷”。稳健回归的菜单程序模块，据有独立性、完备性和良好的移植性。

Description

一种稳健遥感单木冠幅与林木直径回归方法

技术领域

本发明是一种稳健遥感单木冠幅与林木直径回归方法，属于计算机程序及测量技术领域。

背景技术

通过遥感影像获取单木冠幅减轻了林木调查的工作量，由单木冠幅与林木直径的回归模型，可以达到估测林木直径的目的。

因为树冠受各种随机影响大，从遥感测量单木冠幅与实地测量林木直径的线性关系中，有许多异常值点的存在。特别是林分密度作用不易控制。而通常采用的最小二乘法，对异常值点较敏感，这是因为最小二乘法是用残差平方和的极小值“优化”估计参数的。当数据中有异常值点，残差也会异常大，它对估计参数的影响呈平方倍数增加，从而使得线性模型的参数估计不稳定。在数理统计中，人们用模型“稳健性”刻画一个模型受异常值点影响的程度。

为解决这个问题，通常采用2-3倍标准差剔除异常值点。它是以小概率事件在一次抽样中不可能出现为前提条件来删除异常值点(异常值点不是错误点，错误点无论在什么条件下均应删除)。存在两种情况，一种情况是用2-3倍标准差剔除异常值点，可能剔不尽异常值点，还可能剔除正常点；另一种情况是观测数据很重要不能删除“异常值点”，例如润扬大桥的锚锭的冻胀力的21个(日)观测数据出现了1个“异常值点”，由于不是人为因素，无法删除它，因此采用Huber提出的M-估计法，也就是稳健估计(Robust Estimatours)的方法，即不剔除异常值点的稳健回归方法。

发明内容

为了克服现有技术的不足,本发明提供一种稳健遥感单木冠幅与林木直径回归方法。

一种稳健遥感单木冠幅与林木直径回归方法,含有以下步骤：

步骤1)、单木冠幅与林木直径相关数据获取；

步骤2)、稳健回归模型的构建；

步骤3)、单木冠幅与林木直径的稳健回归模型参数的迭代确定。

其中步骤1)冠幅与林木直径相关数据获取有以下两种方法：

方法A)、遥感影像冠幅提取法：

在野外选一块林木中等郁闭、中等大小的林地作为试验场地，用差分GPS测定林内“相关树”林木的地理坐标位置，然后测量相应位置的林木直径。将GPS测定的地理坐标编成一个文件，一次性加载到Arcgis系统，与高分辨率遥感影像叠加，再用遥感影像的单木树冠提取技术测量它的冠幅，从而形成遥感影像冠幅与林木直径的相关数据；

方法B)、目视识别法：

在林地边缘找一个明显地物标，打引线进入林内样地，首先识别样地的角点，再从遥感影像中找到样地,通过遥感影像目视识别法，勾绘样地内“相关树”的林木树冠轮廓，测量冠幅,并记录已测冠幅的林木地理坐标。根据记录的林木地理坐标，在地面样地测量相应的林木直径，也可获取遥感影像冠幅与林木直径的相关数据。

其中步骤2)稳健回归模型的构建有以下步骤：

对于多元线性回归模型：Y＝Xβ+e,

式中：Y是n×1的因变量观测值向量，这里的因变量为林木直径的测量值(单位是cm)向量；

X是n×P的自变量观测值矩阵，P为维数，这里的自变量为单木冠幅遥感测量值(单位是m)的一维矩阵；

β是P×1的待估参数向量；

e是n×1的误差向量；

确定参数β通常用最小二乘法(LS),也就是使估计残差平方和Q_LS(称为目标函数)为最小：

式中：y_i是第i个样本的因变量观测值；

x⁽ⁱ⁾是第i个样本的，P维自变量观测向量；

是第i个样本的残差的平方；

从而获得该线性模型的最优无偏估计：

式中：y是n×1的因变量观测值向量；

x是n×P的自变量观测值矩阵，x'是它的转置矩阵，(x'x)^-1是可逆矩阵x'x的逆矩阵；

是P×1的最小二乘估计参数向量；

最小二乘法的优点是使线性模型的解有良好的解析性质,它的缺点是模型受异常值的影响进行了平方的放大，引起解的不够稳定。

为了减小异常值对线性模型的影响程度，采用Huber的M-估计法，构造稳健回归模型，也就是对原目标函数引入新的权重函数w_i(如附图3):

其中，k为调和常数，这里取k＝1.345,

代入，则新的目标函数为：

它们对残差起到抑制作用，由此获得第一次稳健回归的解：

同理可获得第s次稳健回归的解：

式中：w^(s)是第s次引入的权重函数，它是s-1次线性回归参数估计残差的函数；用迭代方法获得稳健回归的最终解。

本发明的优点：

在典型的东北中温带针阔叶混交林区，在高空间分辨率的单木树冠影像提取基础上，为克服传统最小二乘法对“异常值点”的数据过于敏感的“缺陷”，采用Huber提出M-估计法，成功地构造了遥感影像单木冠幅与林木直径稳健回归模型。

根据M-估计法引进目标函数的加权因子分布，从而改变传统的最小二乘法各样本(包括异常点)权重等同的“缺陷”。M-估计法中权重因子是随着样本的残差而变化，残差越大的，权重就越小，反之亦然。这样就削弱了残差平方和的放大作用。根据稳健回归的散点图，稳健回归抑制回归直线向上漂移的趋势。

稳健回归的菜单程序模块，据有独立性、完备性和良好的移植性。

附图说明

当结合附图考虑时，通过参照下面的详细描述，能够更完整地理解本发明，但此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定，如图其中：

图1为本发明的稳健回归模型构建的流程图。

图2为本发明的迭代计算稳健回归估计值流程图。

图3为本发明的权重与残差关系图。

下面结合附图和实施例对本发明进一步说明。

具体实施方式

显然，本领域技术人员基于本发明的宗旨所做的许多修改和变化属于本发明的保护范围。

实施例1：如图1、图2、图3所示，一种稳健遥感单木冠幅与林木直径回归方法,一般分为以下几步：第1步单木冠幅与林木直径相关数据获取；第2步稳健回归模型的构建；第3步单木冠幅与林木直径的稳健回归模型参数的迭代确定；

具体含有以下步骤：

步骤1)单木冠幅与林木直径相关数据获取采用以下方法：

方法A)、遥感影像冠幅提取法：

在野外选一块林木中等郁闭、中等大小的林地作为试验场地，用差分GPS测定林内“相关树”林木的地理坐标位置，然后测量相应位置的林木直径。将GPS测定的地理坐标编成一个文件，一次性加载到Arcgis系统，与高分辨率遥感图像叠加，再用遥感影像单木树冠提取技术测量它的冠幅，从而形成遥感影像冠幅与林木直径的相关数据；

方法B)、目视识别法：

在林地边缘找一个明显地物标，打引线进入林内样地，首先识别样地的角点，再从遥感影像中找到样地,通过遥感影像目视识别法，勾绘样地内“相关树”的林木树冠轮廓，测量冠幅,并记录已测冠幅的林木地理坐标。根据记录的林木地理坐标，在地面样地测量相应的林木直径，也可获取遥感影像冠幅与林木直径的相关数据。

步骤2)稳健回归模型的构建有以下步骤：

由单木冠幅与林木直径的测量资料(x_i,y_i),i＝1,2,......n建立的多元线性回归模型为：

Y＝Xβ+e

式中：Y是n×1的因变量观测值向量，这里的因变量为林木直径的测量值(单位是cm)向量；

X是n×P的自变量观测值矩阵，P为维数，这里的自变量为单木冠幅遥感测量值(单位是m)的一维矩阵；

β是P×1的待估参数向量；

e是n×1的误差向量；

确定参数β的最小二乘(LS)解是使估计残差平方和Q_LS(称为目标函数)为最小：

式中：y_i是第i个样本的因变量观测值；

x⁽ⁱ⁾是第i个样本的，P维自变量观测向量；

是第i个样本的残差的平方，显然它对Q_LS起着平方放大作用；

β是P×1的估计P维参数向量；

从而获得该线性模型的解析解即最优无偏估计：

式中：y是n×1的因变量观测值向量；

x是n×P的自变量观测值矩阵，x'是它的转置矩阵，(x'x)^-1是可逆矩阵x'x的逆矩阵；

是P×1的最小二乘估计参数向量；

为了减小异常值对模型的影响程度，采用Huber提出M-估计法构建稳健回归模型，就是对原目标函数引入新的权重函数w_i，得到新的目标函数Q_new：

式中：w_i是第i个样本的权重；

是第i个样本的经过加权重的残差平方和；

当权重w_i≡1时，上式就是通用的最小二乘估计的目标函数。为了改进目标函数，使之少受观测异常值的影响，加权函数w_i须满足以下条件(如图3所示)：

①权重最大值为即当而且

②对称性w_i(e)＝w_i(-e)；

③w_i将随着残差绝对值的加大而减小，即单调性；

满足这3个条件的权重函数有多个，在本步骤中令：

其中，k为调和常数，这里取k＝1.345,

将权重函数w_i代入，则新的目标函数为：

由此获得一次稳健回归的解：

式中：y是n×1的因变量观测值向量；

x是n×P的自变量观测值矩阵；

w⁽¹⁾是第1次引入的权重函数；它是上一次线性回归参数估计残差的函数；

同理可获得第s次稳健回归的解：

式中：w^(s)是第s次引入的权重函数，它是s-1次线性回归参数估计残差的函数；用迭代方法获得稳健回归的最终解。

步骤3)单木冠幅与林木直径的稳健回归模型参数的迭代确定有以下步骤：

步骤A)、首先用常规的最小二乘法估计线性模型的初始参数

式中：w⁽⁰⁾≡1是单位矩阵；

步骤B)、计算第S次迭代的权重函数；

在迭代的第s步，记上一次的迭代的残差e^(s-1)，计算第s步权重并将其按对角阵排列：

显然，权重函数是残差的反比函数，

步骤C)、计算按第s步的权重矩阵W^(s)加权最小二乘法的参数的估计值；

从而形成一次稳健回归的解；

步骤D)、判断(其中，δ为迭代误差，在本步骤中设为0.001)转入第B步，进入下一次迭代；否则返回第E步，跳出迭代；

步骤E)、根据β^(s)构建单木冠幅与林木直径稳健的线性模型。

实施例2：如图1、图2、图3所示，一种稳健遥感单木冠幅与林木直径回归方法，还包括如下步聚：

步聚1)、冠幅资料获取；

由于遥感图像上单株树冠与林木树冠存在映射关系，但不一定能做到一一映射关系，所以只能从地面森林样地中获取存在对应关系的树木(相关树)。按照遥感影像的单木树冠定义，这种相关树只能在林地内选取。

在野外，用差分GPS测定林内林木的位置并测量相应的林木直径等，将GPS测定的地理坐标编成一个文件，可一次性加载到Arcgis系统，与高分辨率遥感图像叠加。再用遥感影像单木树冠提取技术测量它的冠幅。这种作法，理论上可做到遥感图像上单株树冠与地面林木(不含下层林木)的一一映射关系。

步聚2)、稳健回归模型的构建同上；

步骤3)单木冠幅与林木直径的稳健回归模型参数的迭代确定同上；

实施例3：如图1、图2、图3所示，一种稳健遥感单木冠幅与林木直径回归方法，还包括如下步聚：

步骤1)冠幅与林木直径相关数据获取；

由于遥感图像上单株树冠与地面林木树冠的存在映射关系，但不一定能做到一一映射关系，所以只能从地面森林样地中选取存在对应关系的树木。

采用目视识别法：

在林地边缘找一个明显地物标，打引线进入林内样地，首先识别样地的角点，再从遥感影像中找到样地,通过遥感影像目视识别法，勾绘样地内“相关树”的林木树冠轮廓，测量冠幅,并记录已测冠幅的林木地理坐标。根据记录的林木地理坐标，在地面样地测量相应的林木直径，也可获取遥感影像冠幅与林木直径的相关数据。此外，也可通过遥感影像获取冠幅与林木直径的相关数据。

步骤2)稳健回归模型的构建同上；

步骤3)单木冠幅与林木直径的稳健回归模型参数迭代，同上。

如上所述，对本发明的实施例进行了详细地说明，但是只要实质上没有脱离本发明的发明点及效果可以有很多的变形，这对本领域的技术人员来说是显而易见的。因此，这样的变形例也全部包含在本发明的保护范围之内。

Claims (1)

1.一种稳健遥感单木冠幅与林木直径回归方法,其特征在于含有以下步骤：

步骤1)、单木冠幅与林木直径相关数据获取；

步骤2)、稳健回归模型的构建；

步骤3)、单木冠幅与林木直径的稳健回归模型参数的迭代确定；

其中步骤1)冠幅与林木直径相关数据获取有以下两种方法：

方法A)、遥感影像冠幅提取法：

在野外选一块林木中等郁闭、中等大小的林地作为试验场地，用差分GPS测定林内“相关树”林木的地理坐标位置，然后测量相应位置的林木直径，将GPS测定的地理坐标编成一个文件，一次性加载到Arcgis系统，与高分辨率遥感影像叠加，再用遥感影像的单木树冠提取技术测量它的冠幅，从而形成遥感影像冠幅与林木直径的相关数据；

方法B)、遥感影像目视识别法：

在林地边缘找一个明显地物标，打引线进入林内样地，首先识别样地的角点，再从遥感影像中找到样地,通过遥感影像目视识别法，勾绘样地内“相关树”的林木树冠轮廓，测量冠幅,并记录已测冠幅的林木地理坐标，根据记录的林木地理坐标，在地面样地测量相应的林木直径，也可获取遥感影像冠幅与林木直径的相关数据；

其中步骤2)稳健回归模型的构建有以下步骤：

对于多元线性回归模型：Y＝Xβ+e,

式中：Y是n×1的因变量观测值向量，这里的因变量为林木直径的测量值(单位是cm)向量；

X是n×P的自变量观测值矩阵，P为维数，这里的自变量为单木冠幅遥感测量值(单位是m)的一维矩阵；

β是P×1的待估参数向量；

e是n×1的误差向量；

确定参数β通常用最小二乘法(LS),也就是使估计残差平方和Q_LS(称为目标函数)为最小：

$Q_{LS}=\underset{i=1}{\overset{n}{\Σ}}{(y_i-x^{\left(i\right)}\mathrm{\β})}^2=\underset{i=1}{\overset{n}{\Σ}}{e}_i^2$

式中：y_i是第i个样本的因变量观测值；

x⁽ⁱ⁾是第i个样本的，P维自变量观测向量；

是第i个样本的残差的平方；

从而获得该线性模型的最优无偏估计：

$\widehat{\mathrm{\β}}={\left(x^{\′}x\right)}^{-1}{x}^{\′}y$

式中：y是n×1的因变量观测值向量；

x是n×P的自变量观测值矩阵，x'是它的转置矩阵，(x'x)^-1是可逆矩阵x'x的逆矩阵；

是P×1的最小二乘估计参数向量；

最小二乘法的优点是使线性模型的解有良好的解析性质；

为了减小异常值对线性模型的影响程度，采用Huber的M-估计法，构造稳健回归模型，也就是对原目标函数引入新的权重函数w_i:

$w_i=\left\{\begin{array}{cc}\frac{1}{2}& \left|e_i\right|\≤k\\ \frac{k\left|e_i\right|-\frac{1}{2}{k}^2}{e^2}& \left|e_i\right|>k\end{array}\right.$

其中，k为调和常数，这里取k＝1.345,

代入，则新的目标函数为：

$Q_{new}=\underset{i=1}{\overset{n}{\Σ}}{w}_i{e}_i^2=\underset{i=1}{\overset{n}{\Σ}}\left\{\begin{array}{cc}\frac{1}{2}{e}_i^2& \left|e_i\right|\≤k\\ k\left|e_i\right|-\frac{1}{2}{k}^2& \left|e_i\right|\≥k\end{array}\right.$

它们对残差起到抑制作用，由此获得第一次稳健回归的解：

${\widehat{\mathrm{\β}}}^{\left(1\right)}={\left(x^{\′}{w}^{\left(1\right)}x\right)}^{-1}{x}^{\′}{w}^{\left(1\right)}y$

同理可获得第s次稳健回归的解：

${\widehat{\mathrm{\β}}}^{\left(s\right)}={\left(x^{\′}{w}^{\left(s\right)}x\right)}^{-1}{x}^{\′}{w}^{\left(s\right)}y$

式中：w^(s)是第s次引入的权重函数，它是s-1次线性回归参数估计残差的函数；用迭代方法获得稳健回归的最终解。