CN105841804A - 一种植被二向性反射一次散射贡献项的获取方法 - Google Patents

Abstract

本发明公开了一种植被二向性反射一次散射贡献项的获取方法，包括以下步骤：通过植被生理参数和聚集指数，获得孔隙率；利用步骤1的孔隙率，定义并计算出立足于植被元素的植被四分量面积比例；求取太阳直射与天空散射光强的比例；通过叶子和土壤的相关参数和上述步骤获得的孔隙率、植被四分量面积比例、太阳直射与天空散射光强的比例，获得植被二向性反射一次散射贡献项。该方法准确度高，计算简便，适用于各种植被类型，把植被冠层非各向同性的起因归结为三个方面：太阳、目标、传感器三者之间的几何关系；多尺度群聚；太阳直射与天空散射光强的比例。

Description

一种植被二向性反射一次散射贡献项的获取方法

技术领域

本发明涉及遥感图像识别技术领域，具体涉及一种植被二向性反射一次散射贡献项的获取方法。

背景技术

遥感影像是地球表面的“相片”，真实地展现了地球表面物体的形状、大小、颜色等信息。要从遥感影像中反演得到地表植物的结构和生理参数，必须要准确描述植物、太阳与传感器之间的关系。

植被对太阳的短波辐射(波长为0.3μm-2.5μm)的反射属于植被二向性反射。现有技术中通常用双向反射率分布函数描述二向性反射行为，定义为：

其中，L为传感器所接收到的来自目标的反射辐射亮度，μ_iF和E分别为太阳直射与天空漫射对目标构成的辐照度，μ_i＝cosθ_i，θ_i为太阳天顶角。θ_v，分别是太阳方位角，观测天顶角和观测方位角，λ为波长。

经过半个多世纪的努力，到目前为止，获取植被二向性反射的方法可以概括为几何光学法、辐射传输法和计算机模拟法等几类。

几何光学法主要针对离散植被，以森林为主要研究对象，把每一个树冠设想为规则的几何体，且假定树冠表面和背景对外来辐射的反射满足朗伯反射要求。

L＝K_cL_c+K_gL_g+K_tL_t+K_zL_z (15)

其中K_c代表光照树冠在传感器视场中所占面积比例，K_g代表光照地表在传感器视场中所占面积比例，K_t代表阴影树冠在传感器视场中所占面积比例，K_z代表阴影地表在传感器视场中所占面积比例，L_c代表光照树冠的朗伯反射辐射亮度值，L_g代表光照地表的朗伯反射辐射亮度值，L_t代表阴影树冠的朗伯反射辐射亮度值，L_z代表阴影地表的朗伯反射辐射亮度值。1997年四尺度几何光学法提出后(Chen，1997)，几何光学法有了新的进展，该文指出植被元素在树冠内的群聚使树冠表面不再满足朗伯反射要求，式(15)存在的前提受到质疑。而且几何光学法在应用于浓密森林植被时，求取四分量面积比例变得十分复杂。

辐射传输法是借助辐射强度(或辐射亮度、辐射通量密度)表达的能量守恒，不存在阴影和热点。植被二向性反射的最显著的特征就是不存在衍射，存在阴影和热点。为了解决辐射传输法无热点的困惑，提出了种种补救方法，如Kussk热点法(Kussk，1993)，但是Kussk热点法本质上是几何光学法。把几何光学法拼入辐射传输法，实质上严重的破坏了理论推导的一致性原则。

计算机模拟法分为蒙特-卡罗模拟法(Monte Carlo method，MC)和真实场景模拟法。

蒙特-卡罗模拟法立足于对一个光子生命过程的追踪，模拟光子能量在一次碰撞过程中散射与吸收间的能量比例(徐希孺，2004)，而不追究散射光强在空间的分布。蒙特-卡罗模拟法不存在热点。

真实场景模拟法有两个基本组成部分，一是植被冠层结构场景的虚拟模拟，二是模拟电磁波与植被元素间的相互作用过程。真实场景模拟法可以产生热点，但计算量过大，随植被元素增加，计算量呈几何级数增长。

由于上述原因，本发明人对现有的获取植被二向性反射一次散射贡献项的方法进行了深入研究，以便设计出一种计算简单准确的方法，其能够为不同植被类型的二向性反射一次散射贡献项提供统一计算方法，可以准确的描述植被二向性反射的方向特征。

发明内容

为了克服上述问题，本发明人进行了锐意研究，提出植被二向性反射可以且必须表达为一次散射贡献项和多次散射贡献项之和，并设计出一种植被二向性反射一次散射贡献项的获取方法，该方法计算简单，结果精确，能够适用于不同植被类型，可以准确描述植被二向性反射的方向特征。通过此方法可以反演获得准确的地表植被的结构和生理参数，比如叶面积指数。

本发明基于几何光学原理，立足于植被元素——叶片，把植被冠层非各向同性的起因归结为三个方面：①太阳、目标、传感器三者之间的几何关系②多尺度群聚③太阳直射与天空散射光强的比例，从而使得本发明不仅适用于离散植物，还适用于连续植物和行播植物，提高了计算准确度，降低了计算复杂度。根据本发明的植被二向性反射一次散射贡献项的获取方法，可以方便地获得等效平均叶面积指数。

附图说明

图1示出植被二向性反射一次散射贡献项的获取方法的步骤图。

图2-a示出行播作物水平剖面图。

图2-b示出行播作物垂直剖面图。

图3示出聚集指数ξ(θ)与太阳天顶角θ的关系图。

图4示出黑杉树孔隙率随角度变化关系(Chen，1997)。

图5示出视线与射线间的相关性对植被四分量的影响。

图6-a-1示出入射方向和观测方向交叠的无因次拦截率l_c的第一种取值范围图。

图6-b-1示出入射方向和观测方向交叠的无因次拦截率l_c的第二种取值范围图。

图6-c-1，6-c-2，6-c-3分别示出入射方向和观测方向交叠的无因次拦截率l_c的第三种取值范围图。

图7示出重叠面积计算示意图。

图8示出主平面l_c中来自入射方向的无因次拦截率l_c，i和l_c中来自观测方向的无因次拦截率l_c，v示意图。

图9示出非主平面l_c中来自入射方向的无因次拦截率l_c，i和l_c中来自观测方向的无因次拦截率l_c，v示意图。

图10-a示出主平面内红色波段的植被二向性反射中一次散射贡献项受天空散射光比例因子β影响的变化图。

图10-b示出主平面内近红外波段的植被二向性反射中一次散射贡献项受天空散射光比例因子β影响的变化图。

图11示出植被二向性反射一次散射贡献项随聚集指数变化关系图。

具体实施方式

下面通过附图和实施例对本发明进一步详细说明。通过这些说明，本发明的特点和优点将变得更为清楚明确。

在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。

在本发明中，所述孔隙率为双向孔隙率，即入射方向上的孔隙率或观测方向上的孔隙率，所述观测方向和所述视向为同一方向。

根据本发明提供的植被二向性反射一次散射贡献项的获取方法，在步骤1中，通过植被生理参数和聚集指数，获得孔隙率，其中，植被生理参数包括植被元素对光子的平均拦截率，孔隙率通过下式(1)获得。

$P\left(\mathrm{\θ}\right)=e^{-\mathrm{\ξ}\left(\mathrm{\θ}\right)\frac{G}{\mathrm{\μ}}{\mathrm{LAI}}_a}---\left(1\right)$

孔隙率定义为当光子穿越植被冠层时，未被植被元素拦截的概率。当植被元素(主要指叶子)在植被冠层中作随机分布时，孔隙率可借用泊松概率分布公式(16)求得。

$P\left(n\right)=\frac{{\mathrm{\η}}^n{e}^{-\mathrm{\η}}}{n!}---\left(16\right)$

n代表单位面积拦截光子的个数，P(n)代表单位面积拦截n个光子概率，η代表平均拦截率。孔隙率即n＝0时的发生概率：P(0)＝e^-η

对植被而言，植被元素对光子的平均拦截率为植被冠层的孔隙率为：

$P=e^{-\mathrm{\η}}=e^{-\frac{G_{i,v}}{{\mathrm{\μ}}_{i,v}}{\mathrm{LAI}}_a}---\left(17\right)$

G_i及G_v分别为太阳光入射方向(i)与视向(v)的G函数值，其中，G函数为在植被层的某一高度处，当叶面积体密度函数取1时向太阳光入射方向的垂直面的平均投影值或观测方向的垂直面的平均投影值，μ_i＝cosθ_i，μ_v＝cosθ_v，LAI_a为平均真实叶面积指数，平均真实叶面积指数为叶的面积指数和非叶部分的面积指数之和，非叶部分指树枝、树干。事实上，对光子产生拦截作用的不仅有叶子，还有枝干、花、穗等，所以此处的LAI_a就不是严格意义下的叶面积指数，具有等效意义。

本发明中所述的孔隙率是在考虑了重叠效应的基础上获得的，具体来说，所述孔隙率是通过泊松公式计算得到的，由于叶子是在空间内随机分布的，其不可避免地会发生重叠的效应，在光子直线传播条件下，重叠效应使得孔隙率增加，如果在计算孔隙率时没有考虑到重叠效应，得到的结果准确性较低，例如，当时，孔隙率P＝e^-1≌0.3678。如果此时叶子均匀分布于空间，那么其孔隙率应为“0”。所以本发明中提供的获得所述孔隙率的方法能够得到更为精确的孔隙率；

本发明中所述的孔隙率是在考虑了群聚效应的基础上获得的，具体来说，所述孔隙率应用泊松公式的前提条件是“随机分布”，而植被树冠或植被元素在空间的分布有可能偏离“随机分布”，比如“群聚”或者均匀分布(在自然状态下，均匀分布很少出现，它的出现大都与人类干预有关)。由于群聚使孔隙率增大，均匀分布使孔隙率变小，因此，当θ≠0时，式(17)被推广为式(1)。所以本发明中提供的获得所述孔隙率的方法能方便而清晰地考虑群聚效应，获得的孔隙率更加精确。

本发明提供的获得所述孔隙率中获得ξ的方法是实测不同θ角条件下的孔隙率P(θ)和LAI_a(包含了所有拦截光子的植被元素的等效LAI)，并代入公式(1)，便可以求得ξ(θ)的变化规律。以行播作物为例，取垂直垄向为主平面，剖面图如图2-a和图2-b所示。y方向取单位距离，

$u_l=\frac{\mathrm{LA}}{a_1\·H\·1}$

$P\left(0\right)=\frac{a_2}{A}(1-e^{-{\mathrm{GLAI}}_1})+e^{-{\mathrm{GLAI}}_1}$

$\mathrm{\ξ}\left(0\right)=-\frac{1}{{\mathrm{GLAI}}_a}\ln P\left(0\right)$

${\mathrm{LAI}}_1=\frac{A}{a_1}{\mathrm{LAI}}_a$

其中，LA代表总叶面积，A代表像元面积，a₁代表群聚效应下，叶子全部集中到面积为a₁的像元内，a₂＝A-a₁，在遥感数据采集，如扫描成像时，像元是传感器对地面景物进行扫描采样的最小单元。

当θ＝θ_c，称θ_c为临界角。三角形面积1和2分别相等，ξ(θ_c)＝1。

当时，a₁，a₂，A同步增长，H不变，θ_c向方向移动时，ξ(θ)值不变，其中，φ代表入射方位角与观测方位角差值的绝对值。ξ(θ)与θ的关系如图3所示。

以离散植被为例，Chen(1997)实测了黑杉树孔隙率，P(θ)与θ的关系如图4所示，其中，Chen(1997)指Chen J M，Leblanc S G.1997.A four-scale bidirectionalreflectance model based on canopy architecture.Geoscience and RemoteSensing.IEEE Transactions on，35(5)：1316-1337。

图4表明当θ＞15°时，所测P(θ)与泊松公式拟合相近；当θ＝7.5°时，P(7.5°)可用分布纽曼拟合，m₂＝24，其中，m₂代表植物集群平均大小的推理值。由此可见，θ_c接近15°。

可以推测当θ＜θ_c时，ξ＜1；当θ≥θ_c时，此时虽然植被元素的尺度群聚仍然存在，但较大尺度的孔隙已消失。此时纽曼分布的m₂≌1，可用泊松分布近似，故ξ≌1。

本发明提供的获得所述双向孔隙率的方法中获得ξ的方法不仅为我们提供了一个求取平均真实叶面积指数(LAI_a)的好方法，而且证明遥感反演所获得的叶面积指数LAI_a为非传统意义下的LAI，为由多角度孔隙率测值所获得的等效LAI的近似值。求取平均真实叶面积指数是通过植被冠层平均高度H和树冠(或植株)平均间距，算得θ_c的值。利用θ＞θ_c的多个孔隙率实测值去反推LAI_a(假定非叶成份的G值已知)，再由LAI_a的值倒推θ＜θ_c时的值例如，用上述方法获得图4对应的黑杉树等效平均叶面积指数为1.87(假定G＝0.5)。而Chen(1997)提供的传统实测单树叶面积指数为4.5(G＝0.5)，折合成像元平均叶面积指数为1.1451。显然，其间的差别体现了树干等非叶成份对光子的拦截作用，这足以证明遥感反演所获得的平均真实叶面积指数LAI_a为非传统意义下的LAI，为由多角度孔隙率测值所获得的等效LAI的近似值。

根据本发明提供的植被二向性反射一次散射贡献项的获取方法，在步骤2中，利用所述孔隙率，获得立足于植被元素的植被四分量面积比例，其中，光照地表在传感器视场中所占的面积比例通过下式(2)获得，光照叶子的地表投影面积在传感器视场中所占的面积比例通过下式(3)获得，阴影地表在传感器视场中所占的面积比例通过下式(4)获得，阴影叶子的地表投影面积在传感器视场中所占的面积比例通过下式(5)获得，

$K_g=e^{-l_i}\·e^{-l_v}\·e^{l_c}---\left(2\right)$

$K_c=1-e^{-l_c}---\left(3\right)$

$K_z=e^{-l_v}(1-e^{-l_i}\·e^{l_c})---\left(4\right)$

$K_t=e^{-l_c}-e^{-l_v}---\left(5\right)$

本发明中所述的立足于植被元素的植被四分量面积比例是在考虑了其与孔隙率及视线与射线间的相关性的基础上获得的，具体来说，立足于植被元素的植被四分量以地表为计算的基准面，以面积比例的形式(无因次拦截率)表达了每一个事件在一次散射中的发生概率。显然，它们的取值与双向孔隙率及视线与射线间的相关性有关，其相关性可用图5表示。

如图5所示，i，v分别代表光线入射方向与传感器视向。两个椭圆形分别代表植被元素(主要是叶子)在地表像元中投影总面积的比例(排除重叠)。根据K_g的定义以及双向孔隙率与无因次拦截率之间的泊松关系，K_g可表达为，

$K_g=e^{-l_i}\·e^{-l_v}\·e^{l_c}---\left(2\right)$

根据植被四分量的定义，可获得式(3)-(5)。

其中，l_c的的取值有三种可能，用图6的a，b，c表示。

如图6-a-1所示，l_c＝0；如图6-b-1所示，l_c取值居中，部分相关；如图6-c-1、6-c-2、6-c-3所示，l_c取最大值，完全相关(最大相关)。在图6-a-1，图6-b-1，6-c-1、6-c-2、6-c-3三种情况下，均满足l_i+l_v-l_c≥0。。

决定于决定入射方向和观测方向交叠的无因次拦截率l_c的因素有三个：

1.入射方向的无因次拦截率l_i与观射方向的无因次拦截率l_v的大小；

2.入射方向的无因次拦截率l_i与观射方向的无因次拦截率l_v的形状；

3.入射方向的无因次拦截率l_i与观射方向的无因次拦截率l_v之间的距离。

l_i与l_v的大小决定于鉴于植被元素在空间呈随机分布，当θ_i，v＝0，l_i与l_v呈圆形，其半径为

当θ_i，v≠0时，在垂直i或v的平面内仍为圆形。但投影到水平面则呈椭圆形，其短轴仍为r₀，长轴用r_i，v，l表示，其取值为

用A_i与A_v来表示椭圆l_i与l_v中心点在x轴(为椭圆长轴的取向)上的位置，则其中H＝H₁+H₂，H₁代表冠层厚度，H₂代表基底(树干)高，μ_l代表植被元素(叶子)体密度函数。

在有因次空间与无因次空间之间进行转换，就是在两个空间中的数值间建立起一一对应的关系。叶子的浓密程度是表征植被几何特性的关键参数之一，可用单位截面积柱体内叶子总面积值来描述。但它是有因次的，有因次量会因为所取计算单位不同而不同，故引入体密度函数μ_l，

${\mathrm{\μ}}_l=\frac{\mathrm{LA}}{D}$

LA代表柱体内总叶面积值，D代表柱体体积，μ_l仍为有因此量，它的因次为|L|^-1，定义平均真实叶面积指数LAI_a＝μ_l·H，则LAI_a为无因次量，它取值唯一。所以一旦μ_l取值确定，则在有因次量H和无因次量LAI_a间就确立了一一对应的关系。可以解释为在无因次空间取值为“1”时，对应有因次空间应有的取值。所以，

H/(1/μ_l)＝Hμ_l＝LAI_a。

当植被冠层离地时(设H₂为其基底高)，则对应无因次空间的取值应为(H₁+H₂)μ_l。

对于如图6-a-1，l_c＝0。

对于如图6-b-1，需要考虑两种情况：1.太阳入射方向和观测方向在同一个平面(此平面为主平面)；2.当太阳入射方向和观测方向不在同一个平面(此平面为非主平面)。

当太阳入射方向和观测方向在同一个平面(此平面为主平面)时，则首先通过联立两个椭圆方程，求得交点的坐标值(即图8中(x，y)值)。

l_c是由i与v方向两个椭圆的部分截面积，可用l_c，i，l_c，v表示构成，如图8所示。

$l_{c,i}=2\·r_{i,0}\·r_{i,l}\·(\frac{1}{2}{t}_0-\frac{1}{4}\sin 2t)$

$l_{c,v}=2\·r_{v,0}\·r_{v,l}\·(\frac{1}{2}{t}_0-\frac{1}{4}\sin 2t)$

$t={\sin }^{-1}\frac{\left|y\right|}{r_{i,v,o}}$

t₀是用弧度值表示的t值，故当交点x位于a_i和a_v之间，则t₀就按上述方法取值；当交点x位于a_i和a_v的一侧，则需要视两椭圆的相对位置和所求得l_c，i和l_c，v位置的具体情况，对其中的一个t₀需要进行(π-t₀)处理，得到新的t₀。

鉴于水平面上的l_c，i和l_c，v分别被放大了cosθ_i和cosθ_v倍，所以

l_c＝l_c，icosθ_i+l_c，vcosθ_v (14)

当太阳入射方向和观测方向不在同一平面内(此平面为非主平面)时，l_i和l_v在地面上的投影如图9所示。同样，l_i和l_v在地面上的投影也是椭圆形，两个椭圆的交叠部分(即图9中的阴影部分)为l_c，l_c又分为l_c，i和l_c，v。为求解l_c，需要确立两个椭圆的方程表达式。为了计算的简洁明了，采用不同的坐标系，其中，l_i位于平面直角坐标系xoy中，l_v位于平面直角坐标系x′oy′中。为太阳入射方向方位角与观测方向方位角的夹角。除了坐标系不同，其他如椭圆中心点、长短轴等的求解均与主平面的求解方法一样。设两个椭圆中心坐标分别为r和r′，长轴分别为a和a′，短轴分别为b和b′，则有，

$\frac{{(x-r)}^2}{a^2}+\frac{y^2}{b^2}=1$

$\frac{{(x^{\′}-r^{\′})}^2}{a^{\′2}}+\frac{y^{\′2}}{b^{\′2}}=1$

设在平面直角坐标系xoy中有一点(x，y)，如图9所示。该点在平面直角坐标系x′oy′中的坐标为(x′，y′)。为点(x，y)和原点的连线与x轴的夹角；为点(x，y)和原点的连线与x′轴的夹角，则有，

进而有，点(x，y)和点(x′，y′)的转化关系如下，

据式(18)-(20)，有，

联立式(19)、(20)、(23)、(24)求解，得到两个椭圆在平面直角坐标系xoy的交点为(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，(x₄，y₄)。以求解l_c，v为例，利用式(23)和(24)，计算交点(x₁，y₁)在平面直角坐标系x′oy′中的坐标为(x₁′，y₁′)。那么在图9中，观测方向在地面投影的椭圆两侧阴影面积的其中之一——S₁的面积为，

$S_1=2a^{\′}{b}^{\′}(\frac{1}{2}{t}_0-\frac{1}{4}\sin 2t)---\left(31\right)$

其中，

当利用交点坐标求解椭圆两侧阴影面积的时候，若是所求的一侧阴影面积没有超过所在椭圆总面积的一半，则t₀就按上述方法取值；当利用交点坐标求解的阴影面积超过所在椭圆面积一半的时候则需要视两椭圆的相对位置和所求得l_c，i和l_c，v位置的具体情况，对其中的一个t₀需要进行(π-t₀)处理，得到新的t₀。

同理，计算出观测方向投影椭圆另一侧的阴影面积S₂。

利用所求得的在平面直角坐标系xoy四个交点(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，(x₄，y₄)可计算出图9中观测方向投影椭圆中画斜线的梯形面积S_t。

观测方向投影椭圆总面积为S_v，则l_c，v为，

$l_{c,v}=\frac{1}{2}(S_v-S_1-S_2-S_t).$

同理，求出l_c，i。

鉴于水平面上的l_c，i和l_c，v分别被放大了cosθ_i和cosθ_v倍，所以

$l_c=l_{c,i}\cos {\mathrm{\θ}}_i+l_{c,v}\cos {\mathrm{\θ}}_v---\left(15\right)$

根据本发明提供的植被二向性反射一次散射贡献项的获取方法，在步骤3中，求取太阳直射与天空散射光强的比例通过下式(6)获得。

$\mathrm{\β}=\frac{E}{{\mathrm{\μ}}_i{F}_i+E}---\left(6\right)$

本发明中所述太阳直射与天空散射光强的比例能准确的表达天空散射光对植被二向性反射一次散射贡献项的影响，具体来说，由于分子散射相位函数的对称性以及分子散射不构成大气散射的主体，大气分子质量稳定少变、太阳光天顶角的变化、波长的增减，都会引起μ_iF_i与E之间的变化。

当投射到大气上界的太阳辐射强度为常数时，经大气作用之后到达地表的辐射能量表现为两部分之和，其中一部分为直射辐射构成的辐射通量，另一部分为天空散射对地表所构成的辐射通量。在此忽略了大气吸收，大气散射辐射中部分向上穿过大气返回太空，这主要由大气分子散射构成。另一部分向下到达地表构成下行辐照度E。由于分子散射相位函数的对称性，以及分子散射不构成大气散射的主体，大气分子质量稳定少变，所以本发明中返回太空的大气散射辐射量稳定少变，进而到达地表的μ_iF_i与E存在互易关系，即大气气溶胶含量增加，E增加，同时μ_iF_i减少；反之，亦然。

太阳天顶角θ_i的变化亦可引起μ_iF_i与E之间的变化。因为θ_i增加，大气散射路径增长，总散射量增加，μ_iF_i减少。所以实测数据表明，θ_i＝-50°时比θ_i＝-40°时热点强度下降，碗边效应增强(chen，1997)。

当λ＞0.5μm，波长的增减亦能引起μ_iF_i与E的变化，。因为太阳辐射可以用6000K的黑体近似，它与λ^-5成正比。而分子散射与λ^-4成正比，气溶胶大颗粒散射与λ^-3＞＞λ^-1成正比。所以，β_nir＞β_red，故近红外波段的热点强度小于红色波段的，但碗边效应却相反。

因此，μ_iF_i与E的相互影响使植被冠层的二向性反射与β值密切相关。绝不能因为假定天空下行辐射亮度在2π空间取常数而断定天空散射光的变化与二向性反射行为无关，或可以简单地以常数予以消除。所以，本发明中所述太阳直射与天空散射光强的比例能准确的表达天空散射光对植被二向性反射一次散射贡献项的影响。

根据本发明提供的植被二向性反射一次散射贡献项的获取方法，在步骤4中，获取光照叶子反射辐射亮度值、光照地表反射辐射亮度值、阴影叶子反射辐射亮度值和阴影地表反射辐射亮度值，其中，鉴于G＝0.5，对平均状况而言，光照叶子反射辐射亮度值通过式(7)获得，光照地表反射辐射亮度值通过式(8)获得，阴影地表反射辐射亮度值通过式(9)获得，阴影叶子反射辐射亮度值通过式(10)获得。

在阴影面，用代表天空下行辐射亮度，某方向天空对地表的辐射通量贡献为2π空间的平均贡献项为(用Sg表示)，由于不相关，并假定各向同性，则令则阴影地表反射辐射亮度值通过式(9)获得，阴影叶子反射辐射亮度值通过式(10)获得。

根据本发明提供的植被二向性反射一次散射贡献项的获取方法，在步骤5中，通过所述光照叶子反射辐射亮度值、光照地表反射辐射亮度值、阴影叶子反射辐射亮度值、阴影地表反射辐射亮度值、植被四分量面积比例、太阳直射与天空散射光强的比例，获得植被二向性反射一次散射贡献项，其中植被二向性反射一次散射贡献项通过式(12)获得。

本发明中所述的植被二向性反射一次散射贡献项的获取方法不仅考虑了光源-目标-传感器三者之间的几何关系对植被二向性反射一次散射贡献项的影响，还考虑了聚集指数和天空散射光对植被二向性反射一次散射贡献项的影响，因此通过本发明提供的方法获得的结果更加准确，更加符合实际情况。

天空散射光对植被二向性反射一次散射贡献项的影响通过下面的例子说明。

令ξ≌1，G＝0.5，有效叶面积指数为1.87，b＝G_iLAI_a，x＝μ_i，K＝0.5867，太阳入射方向和观测方向在同一平面内(即主平面)，太阳入射的天顶角为-40°，观测天顶角从-75°变到75°，观测角度间隔为15°；红色波段β的值由0.1变到0.7，变化间隔为0.05；近红外波段β的值由0.19变到0.42，变化间隔同样为0.01。通过求解K_c、K_g、K_t、K_z、K_nc和K′_nc，根据公式(12)和表1中设定的六分量反射率，β的变化对植被二向性反射一次散射贡献项的影响如图10-a和10-b所示，其中，K_nc代表光照树干的地表投影面积在视场中所占的面积比例，K′_nc代表阴影树干的地表投影面积在视场中所占的面积比例。

表1

从式(12)来论证这些现象的由来。(1-β)K_gρ_g代表光照地表的反辐射亮度，(1-β)K_cρ_v代表光照叶子的反射辐射亮度，S′βK_zρ_g代表阴影地表的反射辐射亮度，(1-S′)βK_tρ_v代表阴影叶子的反射辐射亮度。当β为小值时，(1-β)就为大值，因此计算得到的光照叶子和光照地表的反射辐射亮度就会大，而计算的阴影叶子的反射辐射亮度和阴影地表的反射辐射亮度就会小；相反，当β为大值时，(1-β)就为小值，计算得到的光照叶子的反射辐射亮度和光照地表的反射辐射亮度就会小，而计算的阴影叶子组分反射率和阴影地表的反射辐射亮度就会大。一般而言，计算得到光照叶子的反射辐射亮度和光照地表的反射辐射亮度会高于阴影叶子的反射辐射亮度和阴影地表的反射辐射亮度；另外，根据椭圆六分量计算方法，太阳入射方向同一侧光照部分占主导，太阳入射方向另一侧阴影部分占主导。综合这几个因素，就能很好地解释上述所产生的现象。

聚集指数对植被二向性反射一次散射贡献项的影响通过下面的例子说明。

以陈(1997)的实验为例，计算了θ_i＝-40°，θ_v＝-40°处值与λ值的关系离散植被的有效叶面积指数LAI_e为1.1451，传统意义下的叶面积指数LAI_a为1.87，树的整体高度为7m，冠层高度为6.5m。聚集指数取值范围为0.1-1，间隔为0.1。利用双向孔隙率求解公式，代入不同的聚集指数计算得到相应的双向孔隙率，进而计算六分量，得到不同聚集指数在热点方向上的植二向性反射一次散射项，为图11所示。显然聚集指数对植被二向性反射一次散射贡献项具有重要意义。

实施例

以下通过范例性实施例进一步描述本发明。

以离散植被黑杉林为例，通过实地测量获得表2，其中H₁代表黑杉林的冠层高，H₂代表黑杉林的基底高，LAI_a代表真实叶面积指数，LAI_e代表有效叶面积指数。，ρ_v，red代表红外波段光照叶子反射率，ρ_g，red代表红外波段光照地表反射率，ρ_v，nir代表近红外波段光照叶子反射率，ρ_g，nir代表近红外波段光照地表反射率，

表2

通过实际测量获得孔隙率P随天顶角变化的关系如表3所示，其中聚集指数ξ通过式(1)计算得到，

表3

当太阳入射天顶角为-40°，红色波段，S′＝0.2153时，植被二向性反射一次散射贡献项ρ¹随观测角的变化如表4所示，

表4

当太阳入射天顶角为-40°，近红外波段，S′＝0.333828时，植被二向性反射一次散射贡献项ρ¹随观测角的变化如表5所示，

表5

当太阳入射天顶角为-55°，红色波段，S′＝0.2153时，植被二向性反射一次散射贡献项ρ¹随观测角的变化如表6所示，

表6

以行播作物玉米冠层为例，通过实测获得表7中的数据，其中a₁代表垄宽，a₂代表垄间宽，H代表冠层高度，代表太阳方位角，代表太阳方位角。

表7

当太阳入射角为-27°，红色波段，S′＝0.0410时，植被二向性反射一次散射贡献项ρ¹随观测角的变化如表8所示，

表8

以上结合了优选的实施方式对本发明进行了说明，不过这些实施方式仅是范例性的，仅起到说明性的作用。在此基础上，可以对本发明进行多种替换和改进，这些均落入本发明的保护范围内。

Claims (10)

1.一种植被二向性反射一次散射贡献项的获取方法，其特征在于，该方法包括以下步骤：

步骤1：通过植被生理参数和聚集指数，获得孔隙率；

步骤2：利用所述孔隙率，获得立足于植被元素的植被四分量面积比例；

步骤3：求取太阳直射与天空散射光强的比例；

步骤4：获取光照叶子反射辐射亮度值、光照地表反射辐射亮度值、阴影叶子反射辐射亮度值和阴影地表反射辐射亮度值；

步骤5：通过所述光照叶子反射辐射亮度值、光照地表反射辐射亮度值、阴影叶子反射辐射亮度值、阴影地表反射辐射亮度值、植被四分量面积比例、太阳直射与天空散射光强的比例，获得植被二向性反射一次散射贡献项。

2.根据权利要求1所述的方法，其特征在于，所述植被生理参数包括植被元素对光子的平均拦截率。

3.根据权利要求2所述的方法，其特征在于，所述孔隙率通过下式(1)获得：

$P\left(\mathrm{\θ}\right)=e^{-\mathrm{\ξ}\left(\mathrm{\θ}\right)\frac{G}{\mathrm{\μ}}{\mathrm{LAI}}_a}---\left(1\right)$

其中，P(θ)代表孔隙率，ξ代表聚集指数，G表示G函数，θ代表太阳天顶角或视向天顶角，μ＝cosθ，LAI_a表示平均真实叶面积指数，平均真实叶面积指数包括叶的面积指数和非叶部分的面积指数，代表植被元素对光子的平均拦截率。

4.根据权利要求1所述的方法，其特征在于，所述植被四分量面积比例包括：光照叶子的地表投影面积在传感器视场中所占的面积比例、光照地表在传感器视场中所占的面积比例、阴影叶子的地表投影面积在传感器视场中所占的面积比例和阴影地表在传感器视场中所占的面积比例。

5.根据权利要求4所述的方法，其特征在于，所述光照地表在传感器视场中所占的面积比例通过下式(2)获得，光照叶子的地表投影面积在传感器视场中所占的面积比例通过下式(3)获得，阴影地表在传感器视场中所占的面积比例通过下式(4)获得，阴影叶子的地表投影面积在传感器视场中所占的面积比例通过下式(5)获得，

$K_g=e^{-l_i}\·e^{-l_v}\·e^{l_c}---\left(2\right)$

$K_C=1-e^{-l_c}---\left(3\right)$

$K_z=e^{-l_v}(1-e^{-l_i}\·e^{l_c})---\left(4\right)$

$K_t=e^{-l_c}-e^{-l_v}---\left(5\right)$

其中，K_g代表光照地表在传感器视场中所占的面积比例，K_c代表光照叶子的地表投影面积在传感器视场中所占的面积比例，K_z代表阴影地表在传感器视场中所占的面积比例，K_t代表阴影叶子的地表投影面积在传感器视场中所占的面积比例，l_i代表入射方向的无因次拦截率，e^-li代表入射方向的孔隙率，l_v代表视向的无因次拦截率，e^-lv代表视向的孔隙率，l_c代表入射方向与视线方向上的无因次拦截率之间的相关性，代表相关系数。

6.根据权利要求1所述的方法，其特征在于，所述求取太阳直射与天空散射光强比例通过下式(6)获得，

$\mathrm{\β}=\frac{E}{{\mathrm{\μ}}_i{F}_i+E}---\left(6\right)$

其中，E代表天空漫射对目标构成的辐照度，μ_iF_i代表太阳直射对目标构成的辐照度，μ_i＝cosθ_i，θ_i为太阳天顶角，F_i代表太阳直射对目标物所构成的辐照度。

7.根据权利要求1所述的方法，其特征在于，所述光照叶子反射辐射亮度值通过下式(7)获得，光照地表反射辐射亮度值通过下式(8)获得，阴影地表反射辐射亮度值通过下式(9)获得，阴影叶子反射辐射亮度值通过下式(10)获得，

L_c＝μ_iF_iρ_v,i-v (7)

L_g＝μ_iF_iρ_v,i-v (8)

L_z＝S′Eρ_g,Ω-v (9)

L_t＝(1-S′)Eρ_v,Ω-v (10)

式中，L_c代表光照叶子反射辐射亮度值，L_g代表光照地表反射辐射亮度值，L_z代表阴影地表反射辐射亮度值，L_t代表阴影叶子反射辐射亮度值，S′代表天空对地表的辐射比例，ρ_g,i-v代表土壤的方向—方向反射率，ρ_v,i-v,代表叶子的方向—方向反射率，ρ_g,Ω-v代表土壤的半球—方向反射率，ρ_v,Ω-v代表叶子的半球—方向反射率。

8.根据权利要求7所述的方法，其特征在于，所述S′通过下式(11)获得，

$S^{\′}={\∫}_0^{\frac{\mathrm{\π}}{2}}{\mathrm{cos\θsin\θe}}^{-\frac{\mathrm{\ξ}\left(\mathrm{\θ}\right){\mathrm{GLAI}}_a}{\cos \mathrm{\θ}}}d\mathrm{\θ}---\left(11\right)$

式中，LAI_a代表真实平均叶面积指数，真实平均叶面积指数包括叶的面积指数和非叶部分的面积指数。

9.根据权利要求1所述的方法，其特征在于，所述植被二向性反射一次散射贡献项是通过下式(12)获得的，

ρ¹＝(1-β)K_gρ_g,i-v+(1-β)K_cρ_v,i-v+S'βK_zρ_g,Ω-v+(1-S')K_tρ_v,Ω-v (12)

式中，ρ¹代表植被二向性反射一次散射贡献项,β代表天空散射光比例因子,S′代表天空对地表的辐射比例,ρ_g,i-v代表土壤的方向—方向反射率,ρ_v,i-v代表叶子的方向—方向反射率,ρ_g,Ω-v代表土壤的半球—方向反射率,ρ_v,Ω-v代表叶子的半球—方向反射率。

10.根据权利要求9所述的方法，其特征在于，叶子的方向—方向反射率ρ_v,i-v等于叶子的半球—方向反射率ρ_v,Ω-v等于叶子的漫反射率ρ_v,L，土壤的方向—方向反射率ρ_g,i-v等于土壤的半球—方向反射率ρ_g,Ω-v等于土壤的漫反射率ρ_g,L。