一种针叶植被冠层反射率计算方法及模型

摘要

本发明公开了一种针叶植被冠层反射率计算方法，包括以下步骤：S1：参数识别；输入模型参数，并将输入的参数初步分为三大类：叶片参数、土壤参数和冠层参数；S2：将步骤S1中的叶片参数输入LIBERTY模型进行单片叶光谱信息模拟，得出针叶叶片的反射率和透射率；S3：根据步骤S1中的土壤参数和冠层参数、步骤S2得到的单个叶片的光谱反射率及透射率，计算消光系数及散射系数，并进一步计算得到SAIL模型的参数；S4：将步骤S3计算得到的参数输入SAIL模型，计算冠层的相关反射因子和反射率；S5：根据S4的反射因子和反射率计算冠层反射率。本发明能适用于连续针叶冠层参数反演。

说明书

技术领域

本发明涉及植被参数反演的方法，具体涉及一种利用物理辐射模型反演针叶植被参数的方法。

背景技术

定量遥感在植被反演中越来越广泛的应用，对高效管理植被生理化参数信息、森林火灾预警、农作物培养等方面具有重要意义。近年来物理反演模型以其稳定性和可移植性强的特点成为了国内外热点研究的植被参数反演模型。其中利用叶片模型和冠层模型的物理耦合模型来反演整个冠层的生化组分含量的方法是个热点思路。叶片物理模型考虑光与叶片相互作用的物理机制及叶片的结构，详细描述了光线在叶片内部的传输过程；冠层物理模型解释了电磁波与冠层、叶片、土壤等的相互作用过程，利用冠层模型可以模拟大量的不同种类的冠层光谱样本，而要通过实验方法获取这样的数据是很难的。将叶片模型耦合到冠层模型中可以反演出冠层的生化组分含量，相较经验或者半经验方法意义更加明确，在模型假设范围内，可以任意时间地点应用，不过存在耗时和依赖模型参数精准性的缺点。

近年来对植被参数反演研究较多的模型中阔叶模型比例偏多，针叶植被反演相对较少，针叶植被参数反演还亟待发展。针叶具有自己特殊的结构:没有明显的栅栏组织，其剖面内几乎都是球形细胞。如何表征和模拟针叶叶片的元素存在着很多问题。而且，针叶叶片的形状和尺寸使得光谱特性测量也很困难，即是在实验室内，也很难测定单叶的反射透过率。LIBERTY模型是针对针叶没有明显栅栏组织、大部分为球形细胞的特点发展起来的，用来模拟针叶簇叶或单片针叶的光谱特性。LIBERTY模型视针叶叶片的细胞为标准圆形细胞，认为针叶是由无数叶细胞堆叠在空气中形成，通过若干层细胞的反射率和透过率的迭代过程求得单个针叶的反射率和透过率。

冠层辐射传输模型的代表模型为SAIL，它考虑了多次散射的因素，假设冠层具有如下性质：冠层水平且无限延伸；冠层组分只考虑叶片；冠层分布是各向同性的。叶片模型的输出，即叶片的反射率和透射率，可以作为SAIL模型的输入参数，再结合土壤反射率、叶面积指数、叶倾角分布、太阳天顶角和方位角，观测天顶角和方位角等相关参数，就可以模拟植被冠层400nm到2500nm的反射率。现有的针叶植被参数反演耦合模型是LIBERTY和5-SCALE模型的耦合，不能适用于连续针叶冠层反演。

因此，有必要设计一种适用于连续针叶冠层反演的方法。

发明内容

本发明所解决的技术问题是，针对反演连续大区域针叶植被参数的问题和简化模型参数和计算过程的问题，本发明提出了一种利用针叶叶片LIBERTY模型与冠层辐射传输模型SAIL模型的耦合模型，该模型将针叶单叶片的反射率与透射率模拟出来并作为针叶冠层模型的输入参数，最后模拟出连续针叶冠层的反射率。

本发明的技术方案为：

一种针叶植被冠层反射率计算方法，包括以下步骤：

S1：参数识别

输入模型参数，并将输入的参数初步分为三大类：叶片参数、土壤参数和冠层参数；其中叶片参数包括针叶结构参数、吸收参数以及生化组分含量；土壤参数包括土壤的光谱反射率；冠层参数包括环境参数和冠层结构参数；

S2：将步骤S1中的叶片参数输入LIBERTY模型进行单片叶光谱信息模拟，得出针叶叶片的反射率和透射率；

S3：根据步骤S1中的土壤参数和冠层参数、步骤S2得到的单个叶片的光谱反射率及透射率，计算消光系数及散射系数，并进一步计算得到SAIL模型的参数；

S4：将步骤S3计算得到的参数输入SAIL模型，计算冠层的相关反射因子和反射率，包括太阳直射方向半球反射率、观测方向半球反射率和双直射反射率；

S5：计算冠层反射率。

所述步骤S2具体包括以下步骤：

S2.1：计算叶片细胞介质总吸收系数k：

k＝d×(C_b+f_hC_h+f_wC_w+f_aC_a+f_lC_l+f_nC_n)

其中，d为平均细胞直径；C_b为基吸收；C_h为叶绿素含量，f_h为叶绿素含量吸收系数；C_w为等效水含量，f_w为水分含量吸收系数；C_a为白化吸收，f_a为白化吸收系数；C_l为木质素加纤维素的含量，f_l为木质素加纤维素的含量吸收系数；C_n为蛋白质含量，f_n为蛋白质吸收系数；

S2.2：计算单个叶片细胞的透射率；

为解出单个叶片细胞的透射率和吸收率以及辐射与叶片细胞的相互作用过程，LIBERTY模型假设叶片内部细胞为球形粒子，细胞球体间是由空气间隔隔开的，细胞的表面根据朗伯余弦定律散射入射光。

(1)假设叶片内部细胞表面为理想漫射体并且服从朗伯余弦定律，反射特性呈现明显的各向异性，即从方向θ出来的反射辐射与cosθ成正比。因此，对于非极化辐射的镜面反射，某个方向θ的反射系数m(θ)由菲涅耳方程确定。外部入射辐射的平均反射系数为m_e，指当光线从低折射指数的介质进入高折射指数的介质，为：

内部入射辐射的平均反射系数m_i，指当光线从高折射指数介质进入低折射指数介质，为：

其中，θ_c为临界角；

(2)假设有单位辐射通量进入细胞，因前面假设了细胞壁为漫射体，那么方向θ上每单位立体角的辐射通量为(1/π)cosθ，单位辐射经过一次内部传播，到达表面的所有辐射M为：

首次从细胞内反射的辐射将为m_iM，而有M-m_iM＝M(1-m_i)的辐射将透出细胞，m_iM第二次在细胞内部传播到达表面的辐射为m_iMM，其中m_iM²m_i被散射，m_iM²(1-m_i)被透射出来，反复多次内部传播到表面以及散射和透出，如图4所示。

则单位辐射中被透出的辐射分量为：

τ′＝M(1-m_i)/(1-m_iM)

而吸收的辐射分量为：

1-τ′＝(1-M)/(1-m_iM)

即单个叶片细胞的透射率为τ′，吸收率为1‐τ′；

S2.3：根据单个叶片细胞的透射率为τ′求出无限厚叶片反射率R；

针叶叶片细胞没有明显的栅栏组织，几乎都是球形细胞。针叶叶片是将细胞紧密的聚合起来，即会得到准无限厚的介质。假设细胞表面是水平的并按层排列，那么除了表层的细胞，所有的细胞都将被反射率为R的物质包围。根据单个叶片细胞的透射率为τ′求出无限厚叶片反射率R。

(1)定义参数x_u，x_a和x_d，分别表示细胞内部辐射向上、向相邻和向下散射的分量，它们将分别被上层的，同层相邻的和下层的细胞截住，如图5所示。通过把侧行辐射x_a分解成上行和下行的辐射的分量，用概率系数x表示从细胞内部出射辐射上行辐射的分量，为：

x＝x_u/(1-x_aτ′)

当吸收系数较低时，kd＜1，并且有x_u+x_d+x_a＝1

(2)入射光一部分在表层细胞壁散射，另一部分进入细胞，进入细胞的辐射中又有一部分从细胞内再次透射出来，光线在细胞内反复进行着这样的过程。经过无限次的相互作用后，最后所有贡献到反射率R的辐射可以分为两部分：一部分(P1)是原始单位辐射经过无限次的反射、再被反射后对R的贡献；另一部分(P2)是经过二次、三次…无限次的反射后又进入表层细胞内部的辐射从表层细胞内下行透过后再被反射再进入再上行透过…如此反复而对反射率R的贡献。

无限厚叶片反射率：

R＝P1+P2

整合上面三式，得：

aR²+bR+c＝0

其中

a＝-m_e-τ′+τ′m_e+xτ′-τ′xm_e

b＝2xm_e²+3τ′xm_e-2xτ′m_e²-2x²τm_e+1

c＝-2xm_e-τ′x+2x²τ′m_e

即可得出无限厚叶片反射率R。

S2.4：迭代计算叶片的反射率和透射率；

叶片的反射率ρ和透过率τ是通过若干层细胞的反射率和透过率的迭代过程求的。

假设叶片是由若干层细胞组成的，每层的光学性质相同。

根据无限厚叶片反射率R的表达式，当单层细胞下面没有其他物质的时候，这时候的反射率：

其中R₁、T₁分别表示单层细胞下面没有其他物质时的反射率和透射率。

当无穷层细胞叠在一起的时候，其无限厚叶片反射率为：

通过上式可以解出T₁为：

已知一层细胞的反射率和透过率，两层细胞的反射率和透过率为:

反复迭代，即可得到任意层细胞的针叶叶片的反射率R_N和透射率T_N；迭代公式为：

令ρ＝R_N；τ＝T_N。

所述步骤S3具体包括以下步骤：

S3.1：求解一般叶倾角分布概率：

SAIL模型最常用的六种叶倾角分布类型有喜平型、喜直型、喜斜型、喜极型、均匀型、球面型。因为针叶的叶倾角分布类型特殊，呈簇叶状，因此SAIL模型中常用的分布类型不适用于针叶叶倾角分布，故采用Campbell椭球分布函数来求解针叶叶倾角密度函数，具体算法如下：

服从椭球叶倾角分布的冠层，叶倾角密度函数等同于椭球叶角密度分布，可用面积表示，并最终得出叶倾角密度函数：

式中，α为叶倾角；g(α)表示叶倾角α的概率密度；χ为椭球水平半轴和垂直半轴的比值(χ＝h/l，h为椭球的水平半轴长，l为椭球的竖直半轴长)。

当χ＝1时，Λ＝2，此时椭球叶倾角密度分布为球面分布；

当χ<1时，有：

ε＝(1-χ²)^1/2

当χ>1时，有：

ε＝(1-χ^-2)^1/2

0.1≤χ≤10的范围可用来描述所有具有单峰性的叶倾角密度分布的冠层；

由此求得平均叶倾角的公式：

一般叶倾角分布概率F(θ)＝∫g(α)dα；S3.2：计算消光及散射系数：

SAIL模型是一个基于K‐M理论的4流9参数线性微分方程组，其微分表达式如下：

其中，x为高度，是微分方程的自变量，这里定义冠层顶高度为0，向上为正方向，在冠层内部x取负值；

要解这个线性微分方程组，首先要求取微分方程的系数，然后才能求解微分方程组。而这些系数可大致分为两类：消光系数和散射系数；

(2)求解消光系数

根据以下公式求解太阳直射方向的消光系数：

根据以下公式求解观测方向的消光系数:

其中，θ_s为太阳天顶角，θ_o为观测天顶角；L′＝lai/h，lau为叶面积指数，h为冠层高度；β_s和β_o为临界角，计算公式为：

当确定计算公式中分母不为0且cosβ的计算结果小于1时直接计算出两个临界角的值；当cosβ的计算结果等于1时，两个临介角均等于π；其中cosβ泛指β_s和β_o的余弦；

(2)求解散射系数

(2.1)根据以下公式计算漫辐射E^‐和E⁺的后向散射系数：

漫辐射E^‐和E⁺的前向散射系数计算公式为：

(2.2)根据以下公式计算太阳直射辐射E_S的后向散射系数：

太阳直射辐射E_S的前向散射系数：

(2.3)根据以下公式计算漫辐射E^‐和E⁺的衰减系数：

(2.4)根据以下公式计算观测方向E₀的后向散射系数：

(2.5)根据以下公式计算观测方向E₀的前向散射系数：

(2.6)根据以下公式计算双向散射系数W：

其中，ρ、τ和θ_l分别表示叶片的反射率、透射率及此时相应的叶倾角；ψ为太阳方向与观测方向之间的相对方位角，即两个方向方位角之差；β₁、β₂和β₃为辅助方位角，其取值方法如下：

S3.3计算SAIL模型的参数；

首先，针对S3.2中离散化得到的每一个叶倾角，计算其对应的一般叶倾角分布概率；

然后，将S3.2中离散化得到的各个叶倾角对应的一般叶倾角分布概率分别与步骤S3.2中得到的消光及散射系数相乘后再相加，得到SAIL模型的参数，计算公式为：

Z＝∑F(θ_l)Z(θ_l)

其中，F(θ_l)为叶倾角θ_l对应的一般叶倾角分布概率，Z(θ_l)指代步骤S3.2中得到的k(θ_l)、K(θ_l)、ω(θ_l)、σ(θ_l)、σ′(θ_l)、s(θ_l)、s′(θ_l)、a(θ_l)、v(θ_l)和u(θ_l)中的任意一个，求得的Z相应地指代SAIL模型中的参数k、k、w、σ、σ′、s、s′、a、v、u中的任意一个；

S4：将步骤S3计算得到的参数输入SAIL模型，计算冠层的相关反射因子和反射率；

所述步骤S4具体包括以下步骤：

S4.1：计算叶面积参数；叶面积参数是评价植被冠层反射率的重要指标，其参数计算公式为：

τ_ss＝e^-klai

τ_oo＝e^-Klai

ρ_dd＝(e^mlai-e^-mlai)/(h₁e^mlai-h₂e^-mlai)

τ_dd＝(h₁-h₂)/(h₁e^mlai-h₂e^-mlai)

其中，K和k分别表示观测方向辐射E₀的消光系数和太阳直射辐射E_s的消光系数；lai表示叶面积指数，为模型输入参数；

S4.2：计算加入热点效应后的参数，计算公式如下：

ρ_sd＝C_S(1-τ_ssτ_dd)-D_Sρ_dd

τ_sd＝D_S(τ_ss-τ_dd)-C_Sτ_ssρ_dd

ρ_do＝C_o(1-τ_ooτ_dd)-D_oρ_dd

τ_do＝D_o(τ_oo-τ_dd)-C_oτ_ooρ_dd

ρ_so＝H_o(1-τ_ssτ_oo)-C_oτ_sdτ_oo-D_oρ_sd

其中，各个中间变量的计算方法为：

h₁＝(a+m)/σ

h₂＝(a-m)/σ＝1/h₁

C_S＝[s′(k-a)-sσ]/(k²-m²)

C_o＝[v(K-a)-uσ]/(K²-m²)

D_s＝[-s(k+a)-s′σ]/(k²-m²)

D_o＝[-u(K+a)-vσ]/(K²-m²)

H_s＝(uC_S+vD_s+w)/(K+k)

H_o＝(sC_o+s′D_o+w)/(K+k)

其中，τ_ssoo为热点效应修正参数；由于太阳辐射的问题，特别是稀疏植被区域的土壤和植被的温差相当大，这与其物理表面和气象特点相关，因此需要一个温度修正的过程，同时也可能提供额外的植被冠层信息。热点效应的修正过程为：

根据2/(k+K)效应进行热点纠正，首先根据给定的热点值计算参数：

给定alf的初始值为alf＝10⁶，

1)若热点值q大于0，则认为是有热点效应影响，此时

其中，θ_s为太阳天顶角、θ_o为观测天顶角、为太阳方位角；

若alf计算结果大于200，则令alf＝200；

若alf的计算结果为0，则按以下公式计算出无阴影纯热点效应影响下的参数τ_ssoo和sumint：

τ_ssoo＝τ_ss

若alf的计算结果不为0，则认为是无热点效应影响，按2)中的方法计算τ_ssoo和sumint；

2)若热点值q为0，则认为是无热点效应影响，利用循环相加的方法计算出无热点效应影响下的参数τ_ssoo和sumint，具体步骤为:

第一步、参数初始化：

fhot＝lai*sqrt(K*k)；

x1＝0；

y1＝0；

f1＝1；

fint＝(‐exp(‐alf))*0.05；

sumint＝0；

i＝1

第二步、通过以下迭代计算τ_ssoo和sumint：

步骤1、判断i的取值：

若i小于20，则令x2＝‐log(1‐i*fint/alf；

若i＝20，令x2＝1；

若i大于20，则结束迭代，令τ_ssoo＝f1；

步骤2、依次进行以下计算：

y2＝‐(K+k)*lai*x2+fhot*(‐exp(‐alf*x2))/alf；

f2＝exp(y2)；

sumint＝sumint+(f2‐f1)*(x2‐x1)/(y2‐y1)；

x1＝x2；

y1＝y2；

f1＝f2；

步骤3，令i＝i+1，并转步骤1；

其中，*表示乘法。

S4.3：计算双向反射率：

双向反射率包括两部分，一部分带热点效应，另一部分未带热点效应，两部分求和为最终的计算结果：

ρ_so2＝ρ_so+w*lai*sumint

S4.4：引入土壤反射作用：

土壤位于整个冠层反射系统的最底部，电磁波穿透冠层后会射到土壤地面，经反射后再次反射回冠层顶部，形成冠层反射率的一部分，计算公式为：

d_n＝1－r_s*ρ_dd

S4.5：计算双半球反射率因子：

S4.6：计算太阳直射方向半球反射率：

ρ_sd2＝ρ_dd+(τ_sd+τ_ss)*r_s*τ_dd/d_n

S4.7：计算观测方向半球反射率：

ρ_do2＝ρ_do+(τ_do+τ_oo)*r_s*τ_dd/d_n

S4.8：计算双直射反射率：

ρ_sod2＝ρ_so+((τ_ss+τ_sd)*τ_do+(τ_sd+τ_ss*r_s*ρ_dd)*τ_oo)*r_s/d_n

ρ_so3＝ρ_so2+ρ_so+τ_ssoo*r_s

其中，r_s是土壤反射率。

所述步骤S5具体包括以下步骤：

首先，根据数据库资料求得太阳直射光线与大气散射光线的基本辐射量，分别记为E_s和E_d；

然后，计算大气散射系数skyl参数：

skyl＝0.847-1.61*sin(90-θ_s)+1.04*sin(90-θ_s)²

最后，计算观测方向测定到的辐射量与入射的辐射量的比值，得到针叶植被冠层反射率，公式为：

一种植被冠层反射率的计算模型，公式如下：

其中，R′为冠层反射率，E_s和E_d分别为太阳直射光线与大气散射光线的基本辐射量，skyl为大气散射系数；ρ_do2和ρ_so3分别为观测方向半球反射率和直射反射率；模型参数根据上述植被冠层反射率的计算方法求解。

有益效果：

本发明提出了一种利用针叶叶片LIBERTY模型与冠层辐射传输模型SAIL模型的耦合模型，该模型将针叶单叶片的反射率与透射率模拟出来并作为针叶冠层模型的输入参数，最后模拟出连续针叶冠层的反射率。本发明在充分利用可获取参数的情况下取消植被冠层反射率模拟过程中叶片反射率及透射率输入的过程，直接以叶片模型模拟算法加入冠层SAIL模型以完成耦合过程，通过建立耦合模型模拟出连续针叶冠层的反射率，有效简化连续针叶冠层的反射率模拟过程上参数获取问题并优化算法，加速计算过程，同时耦合模型有利于植被参数反演。

附图说明

图1为本发明原理图；

图2为本发明输入的模型参数；

图3为本发明耦合主体流程；

图4为本发明中辐射传播原理；

图5为本发明细胞内部辐射向上、向相邻和向下散射的分量；

图6为本发明实验模拟的反射率和透射率与实测的反射率和透射率的对比。

具体实施方式

以下结合附图对本发明进行进一步说明。

本发明提出了一种针叶叶片模型与冠层辐射传输模型SAIL模型的耦合模型。该模型首先将给定的参数进行分类，分为叶片参数、冠层参数和土壤参数三部分，之后通过分析计算单个针叶细胞的透射率和吸收率，然后利用分层理论导出给定结构参数下针叶叶片的光谱信息；同时，对于冠层参数首先利用给定的叶倾角分布模型参数计算叶型分布要素，再结合计算出的几何要素计算冠层模型的消光系数及散射系数并与叶型分布模型结合，经过循环过程得到完整叶片分布下的成套冠层参数，最后，利用求得的两个辐射参数及天空散射系数求定冠层的反射率。

本实施案例以南方丘陵地带某颗生长情况良好且冠层光谱容易测定的针叶松树为例对本发明的技术方案进行进一步说明。如附图所示，本发明所涉及的针叶植被冠层反射率的计算方法如下：

经过测定，该树的参数如下表：

叶绿素含量200胡萝卜素含量0等效水厚度0.01叶面积指数4.30干物质含量0.04太阳天顶角30.00叶片结构参数0观测天顶角10.00叶片厚度1.6相对方位角0.00细胞平均直径40热点8.00土壤反射率干土壤基吸收0.045

选择实验波长为400nm-2399nm

S1：参数识别与分类：在得到植被的相关参数后分为三部分，其中叶片参数包括叶片的色素含量、等效水厚度、叶片厚度、细胞平均直径、干物质含量及最大入射角(一般默认为59度)；冠层参数包括叶倾角分布参数、叶面积指数、太阳天顶角、观测天顶角、相对方位角及热点；土壤的背景反射率输入模型参数，

S2：将步骤S1中的叶片参数输入LIBERTY模型进行单片叶光谱信息模拟，得出针叶叶片的反射率和透射率；

(1)计算针叶叶片细胞介质总吸收系数k：

k＝d×(C_b+f_hC_h+f_wC_w+f_aC_a+f_lC_l+f_nC_n)

(2)计算单个叶片细胞的透射率：

为解出单个叶片细胞的透射率和吸收率以及辐射与叶片细胞的相互作用过程，LIBERTY模型假设叶片内部细胞为球形粒子，细胞球体间是由空气间隔隔开的，细胞的表面根据朗伯余弦定律散射入射光。

(3)假设叶片内部细胞表面为理想漫射体并且服从朗伯余弦定律，反射特性呈现明显的各向异性，即从方向θ出来的反射辐射与cosθ成正比。因此，对于非极化辐射的镜面反射，某个方向θ的反射系数m(θ)由菲涅耳方程确定。外部入射辐射的平均反射系数为m_e，指当光线从低折射指数的介质进入高折射指数的介质，为：

内部入射辐射的平均反射系数m_i，指当光线从高折射指数介质进入低折射指数介质，为：

其中，θ_c为临界角；

(4)假设有单位辐射通量进入细胞，因前面假设了细胞壁为漫射体，那么方向θ上每单位立体角的辐射通量为(1/π)cosθ，单位辐射经过一次内部传播，到达表面的所有辐射M为：

首次从细胞内反射的辐射将为m_iM，而有M-m_iM＝M(1-m_i)的辐射将透出细胞，m_iM第二次在细胞内部传播到达表面的辐射为m_iMM，其中m_iM²m_i被散射，m_iM²(1-m_i)被透射出来，反复多次内部传播到表面以及散射和透出，如图4所示。

则单位辐射中被透出的辐射分量为：

τ′＝M(1-m_i)/(1-m_iM)

而吸收的辐射分量为：

1-τ′＝(1-M)/(1-m_iM)

即单个叶片细胞的透射率为τ′，吸收率为1‐τ′；

(3)根据单个叶片细胞的透射率为τ′求出无限厚叶片反射率R；

针叶叶片细胞没有明显的栅栏组织，几乎都是球形细胞。针叶叶片是将细胞紧密的聚合起来，即会得到准无限厚的介质。假设细胞表面是水平的并按层排列，那么除了表层的细胞，所有的细胞都将被反射率为R的物质包围。根据单个叶片细胞的透射率为τ′求出无限厚叶片反射率R。

(3)定义参数x_u，x_a和x_d，分别表示细胞内部辐射向上、向相邻和向下散射的分量，它们将分别被上层的，同层相邻的和下层的细胞截住，如图5所示。通过把侧行辐射x_a分解成上行和下行的辐射的分量，用概率系数x表示从细胞内部出射辐射上行辐射的分量，为：

x＝x_u/(1-x_aτ′)

当吸收系数较低时，kd＜1，并且有x_u+x_d+x_a＝1

(4)入射光一部分在表层细胞壁散射，另一部分进入细胞，进入细胞的辐射中又有一部分从细胞内再次透射出来，光线在细胞内反复进行着这样的过程。经过无限次的相互作用后，最后所有贡献到反射率R的辐射可以分为两部分：一部分(P1)是原始单位辐射经过无限次的反射、再被反射后对R的贡献；另一部分(P2)是经过二次、三次…无限次的反射后又进入表层细胞内部的辐射从表层细胞内下行透过后再被反射再进入再上行透过…如此反复而对反射率R的贡献。

无限厚叶片反射率：

R＝P1+P2

整合上面三式，得：

aR²+bR+c＝0

其中

a＝-m_e-τ′+τ′m_e+xτ′-τ′xm_e

b＝2xm_e²+3τ′xm_e-2xτ′m_e²-2x²τm_e+1

c＝-2xm_e-τ′x+2x²τ′m_e

即可得出无限厚叶片反射率R。

(4)迭代计算叶片的反射率和透射率；

叶片的反射率ρ和透过率τ是通过若干层细胞的反射率和透过率的迭代过程求的。

假设叶片是由若干层细胞组成的，每层的光学性质相同。

根据无限厚叶片反射率R的表达式，当单层细胞下面没有其他物质的时候，这时候的反射率：

其中R₁、T₁分别表示单层细胞下面没有其他物质时的反射率和透射率。

当无穷层细胞叠在一起的时候，其无限厚叶片反射率为：

通过上式可以解出T₁为：

已知一层细胞的反射率和透过率，两层细胞的反射率和透过率为:

反复迭代，即可得到任意层细胞的针叶叶片的反射率R_N和透射率T_N；迭代公式为：

令ρ＝R_N；τ＝T_N。

S3：计算叶倾角分布、消光系数及散射系数

(1)求解一般叶倾角分布概率：

针叶叶倾角分布服从椭球叶倾角分布的冠层，叶倾角密度函数等同于椭球叶角密度分布，可用面积表示，并最终得出叶倾角密度函数：

由此求得平均叶倾角的公式：

(2)几何要素计算：

(3)消光及散射系数计算：

阴影补偿：

①与消光和散射相关的几何因子计算：首先给叶倾角分布加权并离散化，得到一组离散化的中心角度值叶倾角分布。之后，每个叶倾角对应计算一组消光与散射因子，单次计算方案如下：

首先求临界角β_s和β_o，计算公式为：

当确定计算公式中分母不为0且cosβ的计算结果小于1时直接计算出两个临界角的值；当cosβ的计算结果大于1时，两个临介角均等于π；其中cosβ泛指β_s和β_o的余弦；

计算出两个临时界角后就可以计算太阳直射方向与观测方向的消光系数，计算公式分别为：

其中，L′＝lai/h，lai为叶面积指数，h为冠层高度；

②用于计算双向散射系数W的辅助方位角β₁、β₂、β₃计算：辅助方位角的计算主要取决于两个临界角的相对不等关系，取值方法如下：

If：β₁β₂β₃ψ≤|β_s-β_o|ψ|β_s-β_o|2π-β_s-β_o|β_s-β_o|＜ψ＜2π-β_s-β_o|β_s-β_o|ψ2π-β_s-β_oψ≥2π-β_s-β_o|β_s-β_o|2π-β_s-β_oψ

其中，ψ为太阳方向与观测方向之间的相对方位角，即两个方向方位角之差；

得到三个辅助方位角后即可计算单叶片反射率与透射率相关的乘子，用于计算双向散射系数W。

③双向散射系数计算：在得到辅助方位角后，配合S2中已计算得到的单叶片反射率和透射率即可计算得到双向散射系数，计算公式为：

其中，ρ、τ、θ_l分别表示叶片的反射率、透射率及此时相应的叶倾角；令S2中的ρ＝R_α(N)；τ＝T_α(N)。

模型参数解算

(4)散射系数的计算(叶片反射率和透射率的加入)：

①漫辐射E‐和E+的后向散射系数计算公式为：

②漫辐射E‐和E+的前向散射系数计算公式为：

③太阳直射辐射E_S的后向散射系数计算公式为：

④太阳直射辐射E_S的前向散射系数计算公式为：

⑤漫辐射E‐和E+的衰减系数计算公式为：

⑥观测方向E₀的后向散射系数计算公式为：

⑦观测方向E₀的前向散射系数计算公式为：

对于任何冠层来说，叶片倾角不会是固定的，是一个从0到90不断增加的连续过程，因此确定整个冠层的模型参数时，计算方法为叶倾角分布函数作为概率函数F(θ)与模型参数乘积后再相加，计算公式为：

Z＝∑F(θ_l)Z(θ_l)

其中，F(θ_l)为叶倾角θ_l对应的一般叶倾角分布概率，Z(θ_l)为单叶片的参数，Z是经过叶倾角修正模型Z＝∑F(θ_l)Z(θ_k)加以修正后的整个冠层的参数；，Z(θ_l)指代步骤1)和步骤2)中得到的k(θ_l)、K(θ_l)、ω(θ_l)、σ(θ_l)、σ′(θ_l)、s(θ_l)、s′-θ_l)、a(θ_l)、v(θ_l)、u(θ_l)中的任意一个；Z相应地指代针对整个冠层的SAIL模型四个微分方程中的所有系数k、K、w、σ、σ′、s、s′、a、v、u中的任意一个。

S4：计算冠层相关反射因子和反射率，包括太阳直射方向半球反射率、观测方向半球反射率和双直射反射率：

(1)叶面积系数的加入：叶面积指数是评价植被冠层反射率的重要指标，其参数计算公式为：

τ_ss＝e^-klai

τ_oo＝e^-Klai

ρ_dd＝(e^mlai-e^-mlai)/(h₁e^mlai-h₂e^-mlai)

τ_dd＝(h₁-h₂)/(h₁e^mlai-h₂e^-mlai)

(2)解的参数计算：加入热点效应后的参数较复杂，计算公式如下：

ρ_sd＝C_s(1-τ_ssτ_dd)-D_sρ_dd

τ_sd＝D_s(τ_ss-τ_dd)-C_Sτ_ssρ_dd

ρ_do＝C_o(1-τ_ooτ_dd)-D_oρ_dd

τ_do＝D_o(τ_oo-τ_dd)-C_oτ_ooρ_dd

ρ_so＝H_o(1-τ_ssτ_oo)-C_oτ_sdτ_oo-D_oρ_sd

其中各个变量的计算方法为：

h₁＝(a+m)/σ

h₂＝(a-m)/σ＝1/h₁

C_S＝[s′(k-a)-sσ]/(k²-m²)

C_o＝[v(K-a)-uσ]/(K²-m²)

D_s＝[-s(k+a)-s′σ]/(k²-m²)

D_o＝[-u(K+a)-vσ]/(K²-m²)

H_s＝(uC_s+vD_s+w)/(K+k)

H_o＝(sC_o+s′D_o+w)/(K+k)

(3)热点效应的修正：由于太阳辐射的问题，特别是稀疏植被区域的土壤和植被的温差相当大，这与其物理表面和气象特点相关，因此需要一个温度修正的过程，同时也可能提供额外的植被冠层信息。热点效应的修正过程为：

根据2/(k+K)效应进行热点纠正，首先根据给定的热点值计算参数：

给定alf的初始值为alf＝10⁶，

1)若热点值q大于0，则认为是有热点效应影响，此时其中，θ_s为太阳天顶角、θ_o为观测天顶角、为太阳方位角；(若α计算结果大于200则其alf值为200)；若alf的计算结果为0则，计算出无阴影纯热点效应影响下的参数τ_ssoo和sumint，计算公式为：

τ_ssoo＝τ_ss

若alf的计算结果不为0，则认为是无热点效应影响，按2)中的方法计算τ_ssoo和sumint；

2)若热点值q为0，则认为是无热点效应影响，利用循环相加的方法计算出无热点效应影响下的参数τ_ssoo和sumint：具体步骤为:

第一步给出循环中的初始参数：

fhot＝lai*sqrt(K*k)；

x1＝0；

y1＝0；

f1＝1；

fint＝(‐exp(‐alf))*0.05；

sumint＝0；

第二步计算参数：

最后令τ_ssoo＝f1；

上述代码中的*表示乘法；

在程序计算中，前19次循环中令中间参数x2根据给定公式计算结果并代入后面的公式中计算出sumint，第20次循环中令中间参数x2值为1并代入后面的公式中计算出sumint，令τ_ssoo等于第20次循环计算出的f1。

(4)双向反射率计算：双向反射率包括两部分，一部分带热点效应，另一部分未带热点效应，两部分求和为最终的计算结果:

ρ_so2＝ρ_so+w*lai*sum

(5)土壤反射作用的引入：土壤位于整个冠层反射系统的最底部，电磁波穿透冠层后会射到土壤地面，经反射后再次反射回冠层顶部，形成冠层反射率的一部分，计算方案为：1－r_s*ρ_dd

(6)双半球反射率因子计算：

(7)太阳直射方向半球反射率计算：

ρ_sd2＝ρ_dd+(τ_sd+τ_ss)*r_s*τ_dd/(1－rs*ρ_dd)

(8)观测方向半球反射率计算：

ρ_do2＝ρ_do+(τ_do+τ_oo)*r_s*τ_dd/(1－rs*ρ_dd)

(9)双直射反射率计算：

ρ_sod2＝ρ_so+((τ_ss+τ_sd)*τ_do+(τ_sd+τ_ss*r_s*ρ_dd)*τ_oo)*r_s/(1－rs*ρ_dd)

ρ_so3＝ρ_so2+ρ_so+τ_ssoo*r_s

S5：计算冠层反射率：

根据数据库资料可以将太阳直射光线与大气散射光线的基本辐射量求得，分别为E_S和E_d，其次计算大气散射系数skyl参数：

sskyl＝0.847-1.61*sin(90-θ_s)+1.04*sin(90-θ_s)²

冠层反射率的计算方法为观测方向测定到的辐射量与入射的辐射量的比值，公式为：

计算出的植被反射率如图6所示。

本次选取的南方丘陵地带植被理化参数及对比所用光谱均为实测值，模型模拟结果利用RMSE确定可用性及精度问题。

经检测，该模型本例模拟的光谱信息均方根误差(RMSE)为0.18956，模拟算法模拟结果可用且精度较高。

所描述的实例是本发明的一部分实例，而不是全部实例，不能理解为对本发明的限制。基于本发明中的实例，本领域普通技术人员在没有做出创新性劳动的前提下所获得的所有其他实施方式都属于本发明的保护范围。