基于哑变量非线性联立方程组模型的林木参数遥感反演

罗朝沁，孙 华 ，林 辉 ，李际平，陈振雄

（1.中南林业科技大学 林业遥感信息工程研究中心， 湖南 长沙 410004；2国家林业局 中南林业调查规划设计院,湖南 长沙 410014）

基于哑变量非线性联立方程组模型的林木参数遥感反演

罗朝沁1，孙 华1，林 辉1，李际平1，陈振雄2

（1.中南林业科技大学 林业遥感信息工程研究中心， 湖南 长沙 410004；2国家林业局 中南林业调查规划设计院,湖南 长沙 410014）

以湖南省攸县黄丰桥林场Worldview-2影像和地面样地调查数据为基础，采用Mean shift算法对影像进行多尺度分割，提取杉木人工林林木冠幅信息，共提取有效林木冠幅227个，并对提取的冠幅边界信息进行平滑处理。分析调查数据中实测冠幅与影像提取冠幅之间的相关性，结合实测胸径、树高与冠幅的关系，应用曲线估计、非线性联立方程组以及基于哑变量的非线性联立方程组分别建立树高和胸径的最优估算模型，并进行了精度评价。结果表明：将树高与胸径作为哑变量，并进行数量化分级建立的影像冠幅与胸径、树高的非线性误差变量联立方程组模型的拟合效果要优于其他2种方法，树高和胸径模型决定系数R2H和R2D分别为0.899和0.913。模型的适用性检验表明，模型的变动系数、平均百分标准误差均在10%以内，具有较强的稳健性。

林业遥感；林木参数遥感反演；非线性联立方程组；Mean shift 分割算法；黄丰桥林场

树高、胸径等林木参数反映林分内部的结构特征，是森林资源调查中最重要的因子，与林分蓄积量、生物量密切相关。树高、胸径等林木参数的获取一般通过两种途径实现，即地面人工调查和遥感定量反演。对于大范围内的林木参数调查，这种常规的获取方式劳动强度大，成本高[1]。利用遥感技术开展大范围内的林木参数反演是一种比较理想的手段。中低分辨率遥感影像进行林木参数的反演一般是通过利用遥感影像波段反射率、植被指数以及其衍生的变量、纹理因子等与地面实测样地信息建立回归模型[2-4]或空间统计学模型来实现林木参数的反演[5-7]，估测结果地域性强，精度不高，不具推广性。

进入21世纪，随着高分辨率遥感影像的推广与应用，从影像上可以看到明显的冠幅信息，对于散生木或者郁闭度较小林分的单株树冠幅可以直接测量[8]，或者依据树影长度、太阳方位角和坡度等因子，用几何光学的方法计算单木树高[9]。对于郁闭度较高的林分，为了从影像上快速识别和分析冠幅特征，并将冠幅信息从遥感影像数据中分离并提取出来，山谷识别法[10-12]、面向对象的分割技术被认为是比较有效的方式。近年来，冠幅胸径模型以及冠幅树高模型的研究受到了国内外的关注，研究表明冠幅与树高、胸径有着良好的数学关系[13-16]，另外还可以通过冠幅对林分生长、林分竞争[17]以及地上生物量进行预测。与此同时，机载激光雷达技术近年来应用于林木参数反演中，得到了较好的应用效果。机载激光雷达可以获得精度较高的树高信息，但对小范围内的林木位置、冠幅信息的提取，其精度与人工调查还有一定的差距，林木胸径只能通过建立树高胸径模型获取。在林分水平上，林分平均高、胸径等参数一般在单木识别的基础上计算推导，或是通过统计拟合的方式选择最优回归模型进行预测[18-20]。

综上所述，国内外对于冠幅的研究取得了较好的研究成果，利用影像提取的冠幅信息与样地实测数据建立数学模型，或利用半方差函数对单木或小范围内的树高、胸径等林木参数进行估测，研究效果较好。但林木参数遥感反演模型以单因子的线性与非线性回归居多，由于实测数据存在测量误差，建立的林木参数反演模型的误差传递比较明显。文中以Worldview2影像为研究对象，采用曲线估计、非线性联立方程组和基于哑变量的非线性联立方程组模型在攸县黄丰桥林场开展树高、胸径的反演研究，并对3种方法所得结果进行对比，旨在找出一种针对高空间分辨率遥感影像林木参数反演的有效方法。

1 材料与方法

1.1 研究区与数据源

项目研究区位于湖南省攸县黄丰桥林场柏市分场，柏市分场介于东经 113°44′～ 113°43′、北纬 27°14′～ 27°24′之间。分场内的树种主要以杉木Cunninghamia lanceolata、松类Pinusspp.为主。研究所用数据源为Worldview-2数据，含有4个分辨率为1.8 m的多光谱波段和1个0.5 m分辨率的全色波段，数据接收时间为2011年8月11日。

1.2 数据处理

1.2.1 样地数据调查

参考国内外有关林木参数遥感反演的研究方法以及样地设置方法，为了便于从影像上识别林木的位置，在Worldview-2影像上选择郁闭度较低的杉木中龄林作为研究对象。样地调查利用Trimble Geo-XH GPS对样地中心点位置进行定位。利用全站仪测量了样地范围内每棵树的位置，采用Trupulse200激光测距仪进行树高测量，并记录每棵树的年龄、胸径、树高、东西和南北两个方向的冠幅、冠高、枝下高信息，同时记录林分的郁闭度。

1.2.2 遥感数据处理

结合林木参数反演的需求，Worldview-2遥感影像处理主要包括遥感图像的辐射定标、大气校正和正射校正等3个部分。辐射定标根据Digital Globe公式发布的辐射定标技术规程进行。采用ENVI软件中FLAASH模块对Worldview-2影像进行大气校正，Worldview-2影像经过大气校正后，影像亮度发生了明显变化，色调较校正前更加丰富，纹理细节信息明显。为了消除因地形起伏而产生的畸变，选用2台Trimble Geo-XH GPS进行控制点信息采集，选取道路交叉口作为控制点，共采集控制点12个，均匀分布在影像范围内。利用GPS采集信息进行差分计算后的结果，结合1∶10 000地形图数字高程模型（DEM）对遥感影像进行正射校正，控制点总体均方差根误差为4.60个像元（2.3 m），满足研究所需的精度。

1.2.3 冠幅信息提取

为了从影像上快速识别和分析冠幅特征，并将冠幅信息从遥感影像数据中分离并提取出来，面向对象的分割技术被认为是一种有效的方式。研究采用均值漂移算法（Mean Shift algorithm）开展树冠信息提取研究。通过控制Mean Shift算法核的带宽参数（hs，hr）有效的将同质区域影像分割到同一对象块中，最后利用最小区域面积（M）参数逐步实现影像的多尺度分割。整个多尺度分割过程采用由林辉、莫登奎和孙华引入的Mean Shift算法开发的遥感影像分析系统软件V2.0进行。经过反复试验和比较，认为hs=10、hr=6、M=20对于研究区杉木冠幅提取具有较好的效果（见图1，图2）。同时对提取的冠幅边界信息进行平滑处理，利用平滑后的227个杉木冠幅与实测冠幅进行匹配，再建立曲线估计模型，从模型建立的情况来看，实测冠幅与遥感影像提取的冠幅存在良好的线性关系（见图3），其决定系数R2=0.646。

图1 Mean Shift分割结果Fig.1 Segmentation result of Mean Shift algorithm

图2 Mean shift冠幅提取结果Fig.2 Extracted results of tree crows using Mean Shift

1.2.4 树高、胸径参数反演

本研究采取近代统计学中联立方程组法构建影像冠幅胸径模型和影像冠幅树高模型。为了表达方便，用D代表实测胸径，H代表实测树高，C代表实测冠幅，CP代表提取的影像冠幅。具体步骤如下：

1）通过地面实测样木数据，构建D-C模型和H-C模型；

图3 平滑后冠幅与实测值拟合效果Fig.3 Values of tree crown width smoothed against estimated fitting values

2）建立基于分割算法提取的林木冠幅与实测冠幅（C-CP）模型；

3）D-C模型和H-C模型分别与C-CP模型联立。

1.2.5 模型精度评价

利用确定系数R2、估计值的标准差（Standard Error of Estimate，ESEE）、变动系数（Coef ficient of Variation，CCV）和平均百分标准误差（Mean Percent Standard Error，EMPSE）4个评价指标对模型进行检验。

式中：yi表示实测值；表示预测值；表示实测值的均值；p为模型参数个数。

2 结果与分析

2.1 曲线估计模型

分别采用实测冠幅和影像提取的平均冠幅建立冠幅胸径和冠幅树高模型。采用曲线估计中的11种模型进行冠幅树高模型、冠幅胸径模型建模，并选出最佳拟合模型，模型拟合结果如表1所示。从表1中可以看出，利用实测冠幅信息建立的冠幅树高、冠幅胸径模型的拟合最佳曲线均为二次方程，两个模型的决定系数R2分别为0.510和0.520。利用影像提取的冠幅信息建立的冠幅树高、冠幅胸径模型的拟合最佳曲线均为幂函数模型，模型的决定系数R2分别为0.340和0.360。

表1 冠幅胸径/树高最佳曲线估计模型Table 1 Optimal curve estimation models for crown-diameter/height

图4 曲线估计模型残差分布Fig.4 Residual distribution of crown-diameter/height models

表2 冠幅胸径/树高模型检验结果Table 2 Examination results of regression models for crown-diameter/height

2.2 非线性联立方程组反演模型

实测冠幅与胸径以及树高数据分析可知，冠幅与两者之间存在较好的二次方程的关系。由图3可知，实测冠幅与影像冠幅存在线性关系，而实测冠幅与胸径、树高之间是二次方程的关系，是一种非线性关系。假定遥感影像冠幅没有度量误差，则可以通过实测冠幅建立影像冠幅与胸径/树高的非线性联立方程组反演模型。非线性联立方程组的严格数学模型和推导方法见文献[21-23]，研究所建模型的非线性联立方程组运算通过ForStat软件得到。从提取的冠幅227个样本中，抽取157个样本进行构建非线性联立方程组反演模型，剩余70个样本用来做模型精度检验，所得冠幅/胸径、冠幅/树高的非线性联立方程组如下：

式中：联立方程组的误差变量指的是冠幅C和树高H。C的残差平方和为33.159，确定系数R2为0.642；H的残差平方和为619.904，确定系数R2为0.355。

式中：联立方程组的误差变量指的是冠幅C和胸径D。C的残差平方和为33.159，确定系数R2为0.642；D的残差平方和为2 112.708，确定系数R2为0.324。

利用样本验证数据来检验非线性联立方程组模型预测值与实际观测值是否存在较好的线性拟合关系，分析残差的分布是否在置信带内。冠幅树高、冠幅胸径模型的检验结果如图5所示。

图5 非线性联立方程组模型残差分布Fig.5 Residual distribution of non-linear simultaneous model

表3 非线性联立方程组模型检验结果Table 3 Examination results of non-linear simultaneous equations

2.3 基于哑变量的非线性联立方程组反演模型

在考虑树高、胸径以及冠幅测量存在人为误差建立的联立方程组没有得到良好的效果时，应该考虑样地内林木胸径以及树高的差异，将树高和胸径视为哑变量，进行数量化分级，建立基于哑变量的非线性联立方程组反演模型。样地实地调查的过程中，样地中林木的胸径、树高差异较大，样地内最高的杉木为18.6 m，最低的仅为7.2 m，最大的胸径为31.6 cm，最小的为10.4 cm，为了便于计算，分别将树高与胸径划分为4个等级，结果如表4所示。

表4 胸径、树高数量化分级Table 4 Quantitative classification of D and H

将树高和胸径分别视为哑变量，在ForStat软件中构建哑变量非线性联立方程组，将模型中所有参数的初始值设置为1，采用牛顿——拉格朗日方法进行计算。建模结果如下：

式中：C表示实测平均冠幅；CP表示影像上提取的冠幅长度；H表示树高。树高的4个数量化等级表示如下：（1）当E1=E2=E3=0时表示树高第一个等级（I）；（2）当E1=1且E2=E3=0时表示第二个等级（II）；（3）当E2=1且E1=E3=0时表示第三个等级（III）；（4）当E3=1且E1=E2=0时表示第四个等级（IV）。模型中，C的残差平方和为33.159，确定系数R2为0.642；H的残差平方和为97.535，确定系数R2为0.899。

式中：C表示实测平均冠幅；CP表示影像上提取的冠幅长度；D表示胸径。胸径的4个数量化等级表示如下：（1）当C1=C2=C3=0时表示胸径的第一个等级（I）；（2）当C1=1且C2=C3=0时表示第二个等级（II）；（3）当C2=1且C1=C3=0时表示第三个等级（III）；（4）当C3=1且C1=C2=0时表示第四个等级（IV）。模型中，C的残差平方和为33.159，确定系数R2为0.642；D的残差平方和为273.141，确定系数R2为0.913。

表5中模型检验结果表明，冠幅树高模型以及冠幅胸径的拟合效果很好，决定系数R2分别为0.885和0.904。而且，模型检验数据估计值的标准差（ESEE）、变动系数（CCV）、平均百分标准误差（EMPSE）的统计结果较非线性联立方程组有很大的改进，误差基本控制在10%以内。图6中残差点的分布随机性较好，绝大部分验证数据的残差分布在置信带内。从上述分析可知，说明基于哑变量的非线性联立方程组建立影像冠幅与树高以及胸径的反演模型效果十分理想。

表5 哑变量非线性联立方程组模型检验结果Table 5 Examination results of non-linear simultaneous equations when dummy variables were used

图6 哑变量非线性联立方程组残差分布Fig.6 Residual distribution of non-linear simultaneous equations when dummy variables were used

3 结论与讨论

现有针对林木参数的遥感反演，多采用多元统计的分析方法，通过建立实测数据与遥感影像提取的变量之间的回归方程来进行林木参数的反演。目前联立方程组模型还没有应用于林木参数反演的研究当中。本研究以攸县黄丰桥林场Worldview2影像和地面样地调查数据为基础，采用Mean shift算法对影像进行多尺度分割提取林木冠幅信息，对提取的冠幅边界信息进行平滑处理，结合实测冠幅、胸径和树高，分别采用曲线估计、非线性联立方程组和基于哑变量的非线性联立方程组模型，建立林木树高和胸径遥感反演模型。对模型预测值与实际观测值是否存在较好的线性拟合关系进行了检验，分析了模型的残差分布，得出以下几点结论。

（1）Worldview2影像多尺度分割平滑后的冠幅与实测冠幅存在良好的线性关系，决定系数R2=0.646，说明Mean shift算法是林木冠幅信息提取的一种有效方式，通过多尺度的分割可以从高分辨率遥感影像获取量化的林分参数。

（2）地面样地实测冠幅（C）、树高（H）、胸径（D）分别建立的C-H模型、C-D模型的拟合最佳曲线均为二次方程，两个模型的决定系数R2分别为0.51和0.52。利用影像提取的冠幅信息（CP）建立的CP-H模型、CP-D模型的拟合最佳曲线均为幂函数模型，模型的决定系数R2分别为0.340和0.360。而考虑树高、胸径以及实测冠幅存在测量误差，建立影像冠幅与实测冠幅、树高以及胸径的联立方程组反演模型的决定系数R2分别为0.355和0.324。与冠幅与胸径、树高直接建立的反演方程相比，并没有得到明显改善。模型检验数据估计值的标准差、变动系数、平均百分标准误差的统计结果不差上下。

（3）分别将树高与胸径进行数量化分级，将树高与胸径作为哑变量建立的影像冠幅胸径/树高非线性误差变量联立方程组模型的拟合效果比较理想，决定系数R2分别为0.899和0.913，模型的变动系数、平均百分标准误差在10%以内。模型的适用性检验表明，基于哑变量的非线性联立方程组模型的残差点分布随机性较好，以0为基准线上下对称分布。说明利用从影像提取冠幅进行胸径和树高的反演是可行的。

（4）非线性联立方程组模型在林分生长模型领域研究较多，将该模型引入到林木参数遥感反演的研究较少，本研究是在非线性联立方程组模型的基础上，对地面实测数据（胸径和树高）进行数量化分级，并视为哑变量与高分辨遥感影像冠幅提取技术相结合进行林木参数的反演，为林木参数遥感反演提供了一种新的思路。

（5）研究所建立的影像冠幅和实测冠幅、树高以及胸径的基于哑变量的非线性联立方程组模型，主要考虑林木冠幅提取的可操作性，通过与曲线估计、非线性联立方程组模型相比，充分证明基于哑变量的非线性联立方程组模型为本次研究的最优模型。

利用基于哑变量的非线性联立方程组模型与高分辨率遥感影像冠幅定量提取相结合开展林木参数反演，比传统基于植被指数及植被指数衍生变量反演模型更为可靠。在以后的研究中，尤其是针对郁闭度较高的林分可以考虑利用空间变异函数和辐射传输的物理模型结合开展林木冠幅信息提取，进一步提高林木参数的估算精度。

[1]谭炳香,李增元,陈尔学,等. Hyperion高光谱数据森林郁闭度定量估测研究[J].北京林业大学学报,2006,28(3)：95-101.

[2]李崇贵,赵宪文.森林郁闭度定量估测遥感比值波段的选择 [J].林业科学,2005,41(4)：72-77.

[3]杜晓明,蔡体久,琚存勇.采用偏最小二乘回归方法估测森林郁闭度[J].应用生态学报,2008,19(2)：273-277.

[4]曾 涛,琚存勇,蔡体久,等.利用变量投影重要性准则筛选郁闭度估测参数[J].北京林业大学学报,2010,32(6)：37-41.

[5]Franklin S, Wulder M, Gerylo G. Texture analysis of IKONOS panchromatic data for Douglas- fir forest age class separability in British Columbia[J]. International Journal of Remote Sensing ,2002, 22(13)：2627-2632.

[6]冯益明,李增元,邓 广.不同密度林分冠幅的遥感定量估计 [J].林业科学 ,2007,43(1)：90-94.

[7]Jon Pasher, Douglas J K. Multivariate forest structure modeling and mapping using high resolution airborne imagery and topographic information[J]. Remote Sensing of Environment,2010, 114： 1718-1732.

[8]林 辉,吕 勇,宁晓波.基于高分辨率卫星图像的立木材积表的编制[J].林业科学,2004,40(4)：33-39.

[9]崔少伟, 范文义, 金 森, 等. 基于树影与快鸟图像的单木树高提取[J]. 东北林业大学学报, 2011, 39(2)： 47-50.

[10]Gougeon F A. A crown following approach to the automatic delineation of individual tree crows in high spatial resolution aerial images[J]. Canadian Journal of Remote Sensing,1995,21(3)： 274-284.

[11]Culvenor D S. TIDA：an algorithm for the delineation of tree crowns in high spatial resolution remotely sensed imagery[J].Computers and Geosciences, 2002, 28： 33-44.

[12]Michael Palace, Michael Keller, Gregory P Asner,et al. Amazon Forest Structure from IKONOS Satellite Data and the Automated Characterization of Forest Canopy Properties[J]. Biotropica,2008, 40(2)： 141-150.

[13]Hann D W. An adjustable predictor of crown profile for stand grown Douglas- fir trees[J]. Forest Science,1999, 45： 217-225.

[14]Gill S J, Biging G S, Murphy E C. Modelling conifer tree crown radius and estimating canopy cover[J]. Forest Ecology and Management, 2000, 126： 405-416.

[15]Marshall D D, Gregory P J, Hann D W. Crown pro file equations for stand-grown western hemlock trees in northwestern Oregon[J].Canadian Journal of Forest Research, 2003, 33： 2059-2066.

[16]Turan Sönmez. Diameter at breast height-crown diameter prediction models forPicea orientalis[J]. African Journal of Agricultural Research, 2009,4(3)：215-219.

[17]Bechtold W A. Crown-diameter predictions models for 87 species of stand-grown trees in the Eastern United States[J]. Southern Journal of Applied Forestry, 2003, 27： 269-278.

[18]Hall S A, Burke I C, Box D O,et al.Estimating stand structure using discrete-return lidar： an example from low density,fire prone ponderosa pine forests[J]. Forest Ecology and Management, 2005, 208： 189-209.

[19]Thomas V, Finch D A, McCaughey J H,et al.Spatial modelling of the fraction of photo synthetically active radiation absorbed by a boreal mixed wood forest using a lidar hyperspectral approach[J]. Agricultural and Forest Meteorology, 2006, 140：287-307.

[20]Jean François Côté, Richard A Fournier, Richard Egli. An architectural model of trees to estimate forest structural attributes using terrestrial LiDAR[J]. Environmental Modelling &Software, 2011, 26： 761-777.

[21]唐守正,李 勇. 生物数学模型的统计学基础[M].北京： 科学出版社,2002.

[22]唐守正,郎奎建,李海奎.统计和生物数学模型计算(ForStat教程)[M].北京：科学出版社,2009.

[23]刘琼阁,彭道黎,涂云燕.基于偏最小二乘回归的森林蓄积量遥感估测[J].中南林业科技大学学报,2014,34(2)：81-84,132.

Remote sensing retrieval by tree parameters based on dummy variable non-linear simultaneous equations model

LUO Chao-qin1, SUN Hua1, LIN Hui1, LI Ji-ping1, CHEN Zhen-xiong2
(1. Research Center of Forestry Remote Sensing and Information Engineering, Central South University of Forestry and Technology,Changsha 410004, Hunan, China; 2. Central South Forest Inventory and Planning Institute of State Forestry Administration, Changsha 410014, Hunan, China)

A state-owned Huang-Feng-Qiao Forest Farm in Youxian County, Hunan Province was chosen as the study area. Based on Worldview-2 remote sensing data and ground sample survey data of the forest farm, and by adopting Mean shift algorithm, the canopy information of plantation forest of Chinese fir in the farm were extracted with multi-scale segmentation method. Totally, 227 canopy breadth information of effective tested tree were extracted and the tree crown width and crown boundary information extracted were smoothed. The correlation between measured crown width and the crown extracted from the images was studied. By taking into account the relationship between diameter at breast height, tree height and crown width measured, applying curve estimation, nonlinear equations and dummy variable non-linear simultaneous equations, the optimal estimation models of tree height and diameter at breast height were respectively established. Furthermore, the precisions of the model estimation were evaluated. The results show that with theDDBHandHTHas dummy variables respectively, they were graded quantitatively, thus creating a non-linear simultaneous equations model. The fitting effect were much better than those of the first two methods, and the determination coef ficients of tree height and diameter at breast height model (R2H,R2D) were 0.899 and 0.913 respectively. The adaptability testing of the models showed that Standard Error of Estimate, Coef ficient of Variation and Mean Percent Standard Error all were less than 10%, with strong robustness.

forestry remote sensing; remote sensing retrieval by tree pararmeters; non-linear simultaneous equations; mean shift based segmentation algorithm; Huang-feng-qiao Forest Farm in Hubei province

S771.8

A

1673-923X(2015)05-0039-07

10.14067/j.cnki.1673-923x.2015.05.007

2014-06-10

“十二五”国家高技术研究发展计划（863计划）课题（2012AA102001）：“数字化森林资源监测关键技术研究”；中国博士后科学基金项目：林分环境条件下的林木冠幅提取及冠形曲线参数化（2014M562147）；湖南省高等学校科学研究项目：高分辨率遥感影像森林结构参数反演研究（11C1313）

罗朝沁，硕士研究生

孙 华，副教授，博士；E-mail：sunhuayiwen@126.com

罗朝沁,孙 华,林 辉,等. 基于哑变量非线性联立方程组模型的林木参数遥感反演[J].中南林业科技大学学报，2015,35(5)： 39-45.

[本文编校：谢荣秀]