CN102269576A - 一种森林覆盖度及有效叶面积指数的主被动协同反演方法 - Google Patents

Abstract

本发明提供一种森林覆盖度及有效叶面积指数的主被动协同反演方法，该方法包括步骤：1)基于根据机载激光雷达得到的样本区内森林参数，通过Li-Strahler几何光学模型得到每个混合像元的光照地表面积百分比K_g；2)根据每个混合像元的K_g值和光学遥感影像上得到的像元真实反射率，根据线性光谱分解模型获得Li-Strahler几何光学模型中的光照地表的反射率G和其它分量的平均反射率X₀；3)将2)中得到的两个反射率G和反射率X₀，通过Li-Strahler几何光学简化模型得到整个影像上每个真实像元的光照地表面积百分比K_g’，然后得到每个像元的冠层覆盖m和有效植物面积指数PAI_e。本发明利用激光雷达与光学遥感数据相结合，减少混合像元的影响，提高了光学遥感数据反演森林参数的精度。

Description

一种森林覆盖度及有效叶面积指数的主被动协同反演方法

技术领域

本发明涉及一种反演方法，具体来说，涉及一种森林覆盖度及有效叶面积指数的主被动协同反演方法。 

背景技术

在遥感参数反演过程中，混合像元会影响地表生物物理参数的反演精度。由于混合像元是由多种不同覆盖地物的组合，因此混合像元首先影响到像元内任一种地物的覆盖反演(Olthof and Fraser，2007)。例如森林地区的森林植被覆盖，传统的提取植被覆盖度遥感方法是采用植被指数法(Dymond et al.，1992)，由于受到混合像元的影响，这种方法不能提供准确的植被覆盖度信息。而且影响每个混合像元的光谱都是一个多种地物的混合光谱，无论是利用经验回归关系，还是使用模型方法来反演叶面积指数，都会面临着混合光谱分解的问题，如果直接从这种混合光谱中反演其中的叶面积指数，一定会有很大的误差(Hu etal.，2004)，可以说混合光谱分解是提高混合像元反演参数精度的重要一步。随着混合像元分解理论的发展，一些混合像元分解模型被开发出来(Ichoku and Karnieli，1996)，如线性模型、概率模型、几何光学模型、随机几何模型和模糊分析模型。其中一些模型，如线性光谱分解模型、几何光学模型等已被用来反演植被覆盖度(Chopping et al.，2006，2008；Franklin and Strahler，1988；Franklin and Turner，1992；王天星等，2008)和叶面积指数(Hall et al.，1995a，1995b；Hu et al.，2004)，另外线性光谱分解模型和几何光学模型也被同时用于森林参数的反演(Zeng et al.，2008)。但是，在这些模型中，依然没有解决混合像元问题，不仅影响了反演过程，还影响着反演结果的验证。 

另外，要提高传统光学遥感反演森林参数的精度，一个重要途径就是结合激光雷达(LiDAR)数据，但是，单纯利用激光雷达数据来进行成本太高，因此协同激光雷达这种主动光学遥感和传统被动光学遥感 反演森林参数是目前的研究热点。国外已经结合激光雷达数据和传统光学遥感数据开展了很多研究工作，从方法上来归纳，可以把它们分为两类： 

(一)基于图像的融合方法，这种方法主要首先要将激光雷达数据和多光谱图像调成同一分辨率再进行融合，以获得更多的森林参数信息。St-Onge and Achaichia(2001)结合了激光雷达数据和航空摄影测量数据来测量冠层高度，他们用由航空摄影测量得到的海拔高度来代替激光雷达数据的第一个回波，然后用激光雷达地面高度减去航空摄影测量的高度就可得到成像冠层高度模式(CHM，canopy height model)，实验结果表明成像CHM和激光雷达的CHM有很好的相关性。Hudak et al.(2001，2002)融合激光雷达和LANDSAT/ETM+数据来评估森林冠层高度，其中用到的融合方法有非空间spatial方法和空间spatial方法(包括Kriging和Cokriging方法)。Michael and David(2003)结合激光雷达数据和Landsat影像对森林高度进行了调查，根据激光雷达提取的树高与Landsat分割图像上灰度值的关系来反演光学图像的树高。Mundt etal.(2006)将高空间分辨率和高光谱分辨率的数据与激光雷达的分类数据进行了融合，提高了植被覆盖分类的精度。 

(二)基于模型的方法，其中模型包括回归模型和物理模型。而目前大多数的融合研究都是利用回归模型，这种方法只是直接建立由激光雷达数据获得的植被结构参数与光学多光谱遥感数据反演的植被参数之间的关系，并没有从物理机制上说明主动和被动光学遥感数据的结合(Hudak et al.，2006；Lefsky et al.，2005；Popescu et al.，2004；Solberg etal.，2004；Thomas et al.，2006)。而利用物理模型结合激光雷达数据和光学遥感数据来协同反演森林参数的研究并不多，Koetz et al.(2006)同时利用辐射传输模型和激光雷达全波形模型这两种模型来反演森林参数。这种方法首先是利用森林生长模型模拟出辐射传输模型和激光雷达全波形模型所需的参数，然后再根据这些参数分别模拟并建立对应的辐射信号查找表和波形信号查找表，最后将实际获得的遥感影像和激光雷达波形数据代入查找表来反演出森林的生物物理参数。孙国清等(2006)从物理机制上说明了结合大光斑激光雷达数据与多角度光学遥感数据反演森林参数的可行性和优点。总的来说，对于用模型结 合激光雷达点云数据与被动光学遥感数据协同反演的研究工作并不多，模型的选择以及协同方法还需要更多的研究。 

由于机载激光雷达(激光雷达)得到的信息不同于地面上点的实测数据，它得到的是面上信息，有利于与光学遥感数据的协同反演，因此，机载激光雷达可以提取高精度的森林结构参数(包括树高、冠幅、森林覆盖度、叶面积指数等)，这些参数可以作为光学遥感反演叶面积指数的先验知识，也可以作为光学遥感反演结果的验证数据，同时机载激光雷达的分类点云数据(森林植被点云和地面点云)也有利于光学遥感数据中混合像元中端元的提取。但是如何融合光学遥感数据和激光雷达数据，协同两者来反演森林覆盖度和有效叶面积指数是摆在众多研究者面前的一道难题。 

发明内容

对此，本发明提供了一种森林覆盖度及有效叶面积指数的主被动协同反演方法。其从光学遥感反演森林参数中广泛存在着的混合像元分解问题入手，利用从机载激光雷达数据中反演得到光照地表百分比，通过线性光谱分解模型得到混合像元的两个端元信息——光照地表和其他分量，然后结合几何光学模型，反演出森林覆盖率，并在此基础上利用经验模型得到森林的有效叶面积指数。 

本发明的森林覆盖度及有效叶面积指数的主被动协同反演方法包括步骤： 

1)基于根据机载激光雷达得到的样本区内森林参数，通过Li-Strahler几何光学模型得到每个混合像元的光照地表面积百分比K_g； 

2)根据每个混合像元的K_g值和光学遥感影像上得到的像元真实反射率，根据线性光谱分解模型获得Li-Strahler几何光学模型中的光照地表的反射率G和其它分量的平均反射率X₀； 

3)将2)中得到的两个反射率G和反射率X₀，通过Li-Strahler几何光学简化模型得到整个影像上每个真实像元的光照地表面积百分比K_g’，然后得到每个像元的冠层覆盖m和有效植物面积指数PAI_e。 

激光雷达(LiDAR)是一种主动遥感技术，而光学遥感技术(本发明采用的是SPOT-5光学影像)则是被动遥感技术，本发明结合了激光雷 达和SPOT-5数据来进行森林覆盖度和有效叶面积指数的反演，是一种主被动协同的反演方法。 

本发明利用激光雷达与光学遥感数据相结合，只要采集小范围的激光雷达数据，就可以结合光学遥感数据来反演大面积的森林冠层覆盖度及有效叶面积指数。 

本发明通过引入机载激光雷达数据来减少混合像元对光学遥感反演叶面积指数的不确定性问题，提高了光学遥感数据反演森林参数的精度。 

附图说明

图1为本发明的协同反演方法的流程图； 

图2为本发明的利用协同反演技术得到的甘肃省黑河流域大野口地区森林冠层覆盖度； 

图3为本发明的利用协同反演技术得到的甘肃省黑河流域大野口地区森林植被有效叶面积指数； 

具体实施方式

下面以利用激光雷达数据(LiDAR)及光学遥感数据(法国SPOT-5卫星数据)来反演甘肃省黑河流域大野口地区的森林覆盖度及有效叶面积指数为例，详述本发明的反演方法。 

参照图1，在对机载激光雷达提取森林参数和并对光学遥感数据进行预处理之后，开始协同反演。 

Li-Strahler几何光学模型(Li and Strahler，1985；1986；1992)被发展用于遥感影像上像元尺度的树冠大小和树密度的反演，目前已经被广泛用于植被结构参数的反演(Franklin and Strahler，1988；Franklin andTurner，1992；Woodcock et al.，1997；Zeng et al.，2008；李小文和王锦地，1995)。Li-Strahler几何光学模型将遥感影像上的地表反射信号S描述成四个分量的面积百分比权重和： 

S＝K_cC+K_gG+K_tT+K_zZ                  (1) 

其中C是光照树冠的反射信号，G是光照地表的反射信号，T是阴影树冠的反射信号，Z是阴影地表的反射信号；K_c是光照树冠的反 射信号所对应的面积百分比、K_g是光照地表的反射信号所对应的面积百分比、K_t是阴影树冠的反射信号所对应的面积百分比、K_z是阴影地表的反射信号所对应的面积百分比。 

1、混合像元的K_g值确定 

首先根据Li-Strahler几何光学模型和机载激光雷达提供的森林参数(树高、冠幅和冠层覆盖等参数)来得到混合像元中光照地表的面积百分比： 

K_g＝exp{-πm[secθ′_i+secθ′_v-O(θ_i，θ_v，φ)]}        (2) 

其中θ_i、θ_v分别为太阳和卫星的天顶角，φ是太阳与卫星之间的相对方位角；θ′_i和θ′_v分别是从斜视转换到直视坐标系统后的太阳和卫星的天顶角，其转换公式均如6所示；O(θ_i，θ_v，φ)为光照阴影和观测阴影的重叠部分，称为重叠函数，可由下式得到： 

O(θ_i，θ_v，φ)＝(t-sin t cos t)(secθ_i′+secθ_v′)/π  (3) 

其中 

$\cos t=\frac{h}{b}\frac{\sqrt{D^2+{\left(\tan {{\mathrm{\θ}}_i}^{\′}\tan {{\mathrm{\θ}}_v}^{\′}\sin \mathrm{\φ}\right)}^2}}{\sec {{\mathrm{\θ}}_i}^{\′}+\sec {{\mathrm{\θ}}_v}^{\′}}---\left(4\right)$

$D=\sqrt{{\tan }^2{{\mathrm{\θ}}_i}^{\′}+{\tan }^2{{\mathrm{\θ}}_v}^{\′}-2\tan {{\mathrm{\θ}}_i}^{\′}\tan {{\mathrm{\θ}}_v}^{\′}\cos \mathrm{\φ}}---\left(5\right)$

$\tan {\mathrm{\θ}}^{\′}=\frac{b}{r}\tan \mathrm{\θ}---\left(6\right)$

公式4和6中的h、b和r分别是混合像元内的平均树高、树冠的平均长半径和平均短半径(平均冠幅半径)，这些变量都可以从机载激光雷达提取的森林参数和实测数据中得到。 

另外公式2中的另一个参数m是联系遥感影像信号与森林参数的一个重要参数(Woodcock et al.，1997)，当树不互相重叠时，它表示的是样地内的平均冠层覆盖大小，是树密度λ和平均冠幅半径r的函数： 

m＝λr²               (7) 

在求解四个分量在像元中的比重时，都是在假设树在样地中是随机分布的，并不互相重叠，这样冠层覆盖与树密度、冠幅大小的关系可用如下公式来表达(Li and Strahler，1985；Woodcock et al.，1997)： 

CC＝1-e^-πm           (8) 

从公式2可知m在很大程度上影响着K_g及光照地表反射率的求解，最终影响了模型在遥感影像上反演森林参数的精度。 

2、得到每一个真实像元的K_g’值 

这里利用了线性光谱分解模型(方法)，线性光谱分解模型描述的是地表反射率光谱是由多个纯像元的光谱按面积百分比组成。它的数学表达式如下： 

${\mathrm{\ρ}}_j=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{ij}}{s}_i+e_j$

$\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{s}_i=1---\left(9\right)$

其中ρ_ij是第i个纯像元、第j个波段的反射率，s_i是第i个纯像元的面积百分比，e_j是第j个波段的误差。 

公式1把遥感信号看做是四个分量的结合体：光照树冠的反射信号，光照地表的反射信号，阴影树冠的反射信号，阴影地表的反射信号；如果把遥感信号简单的看作是两个分量的结合体，即光照地表的反射信号以及其它物体的反射信号，那么，针对研究内容和研究区的情况，Li-Strahler几何光学模型对混合像元的遥感信号S描述可以简化为一个二分量模型，即公式1可以简化为如下公式： 

S＝K_gG+(1-K_g)X₀          (10) 

其中X₀是混合像元内除光照地表以外的其它总分量，其它参数含义与公式1中的参数相同。 

通过上面计算出的混合像元的K_g值，根据SPOT5影像上得到的像元真实反射率，由线性光谱分解模型，即公式9获得混合像元的光照地表反射率G和其它分量的平均反射率X₀。 

线性光谱分解冲最困难也是最重要的一步就是端元光谱选取，只有得到了端元光谱G和X₀，才能结合大范围的光学遥感数据，利用简化的Li-Strahler模型即公式10计算得到K_g，其原理实际上还是简化Li-Strahler模型。本发明采用因子分析法自动从SPOT5影像上获取端元光谱，因子分析法计算对应端元光谱的抽象矩阵和端元百分比，然后通过建立转换矩阵获取实际的端元光谱及其百分比。 

因子分析模型描述如下： 

1)X＝(x₁，x₂，...，x_p)是可观测随机向量，均值向量E(X)＝0，协方差阵Cov(X)与相关矩阵R相等(只要将变量标准化即可实现)。 

2)F＝(F₁，F₂，...，F_m)(m＜p)是不可测的向量，其均值向量E(F)＝0，协方差矩阵Cov(F)＝I，即向量的各分量是相互独立的。 

3)e＝(e₁，e₂，...，e_p)与F相互独立，且E(e)＝0，e的协方差阵∑是对角阵，即各分量e之间是相互独立的，则模型： 

x₁＝a₁₁F₁+a₁₂F₂+...+a_1mF_m+e₁

X₂＝a₂₁F₁+a₂₂F₂+...+a_2mF_m+e₂

......... 

x_p＝a_p1F₁+a_p2F₂+...+a_pmF_m+e_p

其称为因子分析模型，由于该模型是针对变量进行的，各因子又是正交的，所以也称为R型正交因子模型。其矩阵形式为：x＝AF+e。 

将SPOT5影像各像元的反射率值作为可观测随机向量X，e设为0，通过将各像元端元光谱百分比作为权重系数即a带入矩阵方程中即可求得作为不可测的向量F，即两个端元光谱G和X₀，它们分别代表端元光谱中的光照地表反射率和其它分量的平均反射率。 

在获得像元内两个端元的值后，将这两个端元反射率和SPOT5影像的像元反射率代入Li-Strahlr几何光学模型的简化公式10中，即可反演出SPOT5影像上每一个真实像元的K_g’值。 

3、Li-Strahler几何光学模型反演冠层覆盖和有效PAI 

在获得每一个真实像元的光照地表面积百分比K_g’之后，由公式2的变换式即公式11求得m值，最后由公式8求得像元的冠层覆盖。 

$m=\frac{-\ln \left(K_g\right)}{(\sec {\mathrm{\θ}}_i+\sec {\mathrm{\θ}}_v)(\mathrm{\π}-t+\cos t\sin t)}---\left(11\right)$

另外，还可以根据利用公式12来反演得到森林植被有效植被面积指数，如图3所示。 

PAI_e≈-2ln(e^-πm)＝2πm                (12) 

图2和图3为根据本发明反演出的甘肃省黑河流域大野口地区的森林覆盖度及有效叶面积指数。 

我们利用LAI-2000叶面积测量仪器实地对样区内的一些样点进行了测量，并用该测量结果来验证本专利的协同反演方法，结果如图所 示，本方法与实际测量的真实值之间的R²达到了0.7471，表明本方法具有较高的反演精度。 

Claims (7)

1.一种森林覆盖度及有效叶面积指数的主被动协同反演方法，包括步骤：

1)基于根据机载激光雷达得到的样本区内森林参数，通过Li-Strahler几何光学模型得到每个混合像元的光照地表面积百分比K_g；

2)根据每个混合像元的K_g值和光学遥感影像上得到的像元真实反射率，根据线性光谱分解模型获得Li-Strahler几何光学模型中的光照地表的反射率G和其它分量的平均反射率X₀；

3)将2)中得到的两个反射率G和反射率X₀，通过Li-Strahler几何光学简化模型得到整个影像上每个真实像元的光照地表面积百分比K_g’，然后得到每个像元的冠层覆盖m和有效植物面积指数PAI_e。

2.根据权利要求1所述的森林覆盖度及有效叶面积指数的主被动协同反演方法，其特征在于，步骤1)中的样本区内森林参数包括树高、冠幅和冠层覆盖。

3.根据权利要求1所述的森林覆盖度及有效叶面积指数的主被动协同反演方法，其特征在于，步骤2)还包括步骤：利用因子分析模型计算对应端元光谱的抽象矩阵和端元百分比，然后通过建立转换矩阵获取端元光谱反射率G和其它分量的平均反射率X₀。

4.根据权利要求3所述的森林覆盖度及有效叶面积指数的主被动协同反演方法，其特征在于，其中因子分解模型的转换矩阵为：x＝AF+e，其中x为可观测随机向量，是影像各像元的反射率值，e为0，可观测随机向量，均值向量E(X)＝0，A为权重系数。

5.根据权利要求1所述的森林覆盖度及有效叶面积指数的主被动协同反演方法，其特征在于，步骤3)中的Li-Strahler几何光学简化模型为：

S＝K_gG+(1-K_g)X₀

其中，S表示地表反射信号，G是光照地表的反射信号。

6.根据权利要求1所述的森林覆盖度及有效叶面积指数的主被动协同反演方法，其特征在于，步骤3)中由如下公式求得像元的冠层覆盖m：

$m=\frac{-\ln \left(K_g\right)}{(\sec {\mathrm{\θ}}_i+\sec {\mathrm{\θ}}_v)(\mathrm{\π}-t+\cos t\sin t)}$

其中， $\cos t=\frac{h}{b}\frac{\sqrt{D^2+{\left(\tan {{\mathrm{\θ}}_i}^{\′}\tan {{\mathrm{\θ}}_v}^{\′}\sin \mathrm{\φ}\right)}^2}}{\sec {{\mathrm{\θ}}_i}^{\′}+\sec {{\mathrm{\θ}}_v}^{\′}},$

$D=\sqrt{{\tan }^2{{\mathrm{\θ}}_i}^{\′}+{\tan }^2{{\mathrm{\θ}}_v}^{\′}-2\tan {{\mathrm{\θ}}_i}^{\′}\tan {{\mathrm{\θ}}_v}^{\′}\cos \mathrm{\φ}},$

$\tan {\mathrm{\θ}}^{\′}=\frac{b}{r}\tan \mathrm{\θ},$

h、b和r分别是混合像元内的平均树高、树冠的平均长半径和平均短半径，θ_i、θ_v分别为太阳和卫星的天顶角，φ是太阳与卫星之间的相对方位角，θ′_i和θ′_v分别是从斜视转换到直视坐标系统后的太阳和卫星的天顶角。

7.根据权利要求1所述的森林覆盖度及有效叶面积指数的主被动协同反演方法，其特征在于，步骤3)中由如下公式求得和有效植物面积指数PAI_e：

PAI_e≈-2ln(e^-πm)＝2πm。

Images