CN106909750B - 一种阔叶植被冠层反射率的计算方法 - Google Patents

一种阔叶植被冠层反射率的计算方法 Download PDF

Info

Publication number
CN106909750B
CN106909750B CN201710141821.0A CN201710141821A CN106909750B CN 106909750 B CN106909750 B CN 106909750B CN 201710141821 A CN201710141821 A CN 201710141821A CN 106909750 B CN106909750 B CN 106909750B
Authority
CN
China
Prior art keywords
calculating
blade
reflectivity
canopy
leaf
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710141821.0A
Other languages
English (en)
Other versions
CN106909750A (zh
Inventor
郭云开
安冠星
刘海洋
蒋明
谢琼
周烽松
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changsha University of Science and Technology
Original Assignee
Changsha University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Changsha University of Science and Technology filed Critical Changsha University of Science and Technology
Priority to CN201710141821.0A priority Critical patent/CN106909750B/zh
Publication of CN106909750A publication Critical patent/CN106909750A/zh
Application granted granted Critical
Publication of CN106909750B publication Critical patent/CN106909750B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/04—Inference or reasoning models
    • G06N5/042—Backward inferencing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • General Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Battery Electrode And Active Subsutance (AREA)

Abstract

本发明公开了一种阔叶植被冠层反射率的计算方法及模型，计算方法包括以下步骤：S1：输入参数将输入的参数识别并分类，分为叶片参数、冠层参数及土壤参数；S2：根据叶片参数计算单个叶片的反射率及透射率；S3：根据冠层参数及S2中的叶片参数求得冠层的消光系数及散射系数；S4：根据所求的冠层消光和散射参数求得冠层的相关反射因子和反射率；S5根据冠层的相关反射因子和反射率求得冠层反射率。本发明将PROSPECT模型与SAIL模型耦合，在充分利用可获取参数的情况下取消植被冠层反射率模拟过程中叶片反射率及透射率输入的过程，本发明有效简化植被冠层光谱信息模拟过程上参数获取问题并优化算法，加速计算过程，同时耦合模型有利于植被参数反演。

Description

一种阔叶植被冠层反射率的计算方法
技术领域
本发明涉及高等级公路路域植被健康状况评价，具体涉及一种阔叶植被冠层反射率的计算方法及模型。
背景技术
植被信息反演一直以来都是定量遥感最具前景的研究领域，为开展植被生态环境监测与评价研究提供了有利基础。国内外学者提出了众多植被参数反演模型，主要可分为统计模型和物理模型两大类，其中物理模型以其优秀的广泛适用性及稳定性备受青睐。植被信息的物理反演模型具有相当强的物理基础，不依赖于植被的具体类型和背景环境的变化。然而物理模型在大范围的植被信息反演中存在精度仍然较差且计算量大导致无法满足实际生产应用需求的问题，主要是由于两个原因：一是物理模型对定量遥感的要求很高，需要通过遥感影像反演出地表反射率，同时需要多个难以获得的精确参数作为物理模型的输入参数；另一方面，当前遥感物理模型大多在冠层尺度上，而冠层光谱受到叶片光谱信息，土壤背景反射率及冠层结构等因素的影响，限制了生化参数反演的精度也加大了反演的难度。
近年来研究和应用最多的模型有适用于阔叶的PROSPECT模型、适用于针叶的LIBERTY模型及冠层SAIL模型。PROSPECT模型将阔叶叶片看作是由N 层含有叶片生化物质的粗糙平板和N-1层空气构成的复合平板模型，通过输入各色素浓度、等效水厚度、白物质浓度及叶片结构参数得到这个复合平板模型的反射率和透射率。LIBERTY模型视针叶叶片的细胞为标准圆形细胞，认为针叶是由无数叶细胞堆叠在空气中形成，通过求得单个细胞的光学特性进一步求得无限细胞的光谱信息。冠层SAIL模型是在辐射传输理论的基础上对辐射传输方程的四流近似，模型通过输入太阳天顶角、太阳方位角、观测天顶角、观测方位角、 LAI、土壤反射率、热点及单叶片的反射率与透射率等参数来模拟冠层400nm到 2500nm的反射率。显而易见的是可以将阔叶与针叶模型的输出参数作为SAIL 模型的输入参数来实现模型的耦合，进一步实现植被冠层信息的模拟。因此，国内外学者先后开展了模型耦合及其应用的研究，当前应用较多的是利用 PROSPECT叶片模型与SAIL冠层反射模型耦合，该耦合模型利用修正后的叶片模型与冠层模型模拟植被冠层反射率，具有较高的反演精度和较快的反演速度。然而，当前所用多为发展较为完善的初期模型耦合而成的PROSAIL模型，在未能考虑植被叶片上下表面不对称性及热点效应等问题的情况下对PROSPECT模型与SAIL模型进行耦合，尽管适应性及精度有一定的保证，但在实际应用过程中仍然暴露出了植被参数反演不足等问题，使得耦合模型在实际应用过程中受到一定程序的限制。
因此有必要设计一种新的阔叶植被冠层反射率的计算方法及模型。
发明内容
本发明所解决的技术问题是，针对物理模型在实际应用过程中受限的问题，也为了解决植被物理模型在应用过程中参数获取及反演精度问题，提出一种阔叶植被冠层反射率的计算方法及模型，将单叶片模型利用所获取参数得到的单叶片反射率与透射率作为冠层模型的输入参数，协同冠层模型的其他输入参数模拟植被冠层反射率，耦合后的植被冠层模型可省略部分复杂参数的获取过程，同时耦合后的模型为一体化模型，更加有利于植被参数的反演。
本发明的技术方案为：
一种阔叶植被冠层反射率的计算方法，包括以下步骤：
S1：参数识别；
输入模型参数；并将输入的参数初步分为三大类：叶片参数、土壤参数和冠层参数；其中叶片参数包括色素含量、干物质含量、等效水厚度及叶片结构参数；土壤参数为土壤的光谱反射率；冠层参数为剩余其他参数；
S2：根据步骤S1中的叶片参数输入PROSPECT模型进行单片叶光谱模拟，计算单个叶片的光谱反射率及透射率；
S3：根据步骤S1中的土壤参数和冠层参数计算消光系数及散射系数；
S4：将步骤S3计算得到的消光系数及散射系数输入SAIL模型，计算冠层的相关反射因子和反射率，包括太阳直射方向半球反射率、观测方向半球反射率和双直射反射率；
S5：计算冠层反射率。
所述步骤S2具体包括以下步骤：
S2.1：计算吸收系数K；吸收系数是叶片各种化学物质含量的线性表达，计算公式为：
其中，c(i)为叶片中组分i的浓度，组分i为色素、水、干物质等化学物质；k(i) 为组分i的特定吸收系数；N为叶片结构参数，表示叶片分层数，N可以采用以下经验方程计算得到：
N＝(0.9*SLA+0.025)/(SLA-0.1)
其中，SLA为比叶面积，指每单位干重的叶面积；
S2.2：计算界面处的透射率：
界面处的透射率包括两个：一个是光从空气进入叶片的透射率，一个是光从叶片进入空气的透射率；
自然光线可看作是非极化光，反射后发生完全极化，透射后折射光线发生部分极化。根据Snell-Descarts定律可知道，两种介质接触面上入射角与折射角及两个折射率之间的关系，再根据Schanda给出的电磁波穿过两种介质的透射率即可计算出光线以立体角α由折射率为1的空气入射折射率为n的叶片的平均透射率tav(α,1,n)，令t12＝tav(α,1,n)；
根据Stern的研究成果，光线从折射率为n的叶片射入折射率为1的空气的透射率t21与t12之间存在关系式t21＝n－2·t12，据此计算出界面处的透射率；界面处的透射率的计算为现有技术，可使用现有软件完成计算；
S2.3：计算光线透过平板介质的总透射率τ1；
对于自然光，其透过平板介质的总透射率τ1为整个半球2π整空间内的积分，与吸收系数K及叶片厚度D相关，计算公式为：
其中，t为中间变量；
S2.4：计算单层平板的反射率Rα(1)及透射率Tα(1)：
进行单层平板光谱模拟，单层平板的光谱是整个叶片光谱模拟的核心，主要包括计算单层平板的反射率及透射率，其计算公式为：
S2.5：计算单片叶反射率Rα(N)及透射率Tα(N)：
进行叶片光谱信息模拟，Stokes对光线穿透有限层平板后的光学现象进行了深入研究，得出了光线穿过N层相同反射率和透射率的平板后的总体反射率和透射率的计算理论，此时应用其理论可得出叶片光谱信息的模拟公式：
其中：
所述步骤S3具体包括以下步骤：
S3.1阴影补偿：
①计算与消光和散射相关的几何因子：首先根据一般叶倾角分布概率加权并离散化，得到一组离散化的叶倾角；之后，计算每个叶倾角对应的消光与散射因子；对于叶倾角θl，计算方法如下：
首先，求临界角βs和βo，计算公式为：
其中，为水平面法线方向和叶片法线方向夹角即叶倾角θs为太阳天顶角、θo为观测天顶角；
当确定计算公式中分母不为0且cosβ的计算结果小于1时直接计算出两个临界角的值；当cosβ的计算结果等于1时，两个临介角均等于π；其中cosβ泛指βs和βo的余弦；
然后，计算叶倾角为θl的单个叶片太阳直射方向与观测方向的消光系数，计算公式分别为：
其中，L′＝lai/h，lai为叶面积指数，h为冠层高度；
②计算辅助方位角β1、β2、β3：辅助方位角的计算主要取决于两个临界角的相对不等关系，取值方法如下：
If： β1 β2 β3
ψ≤|β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 ψ
其中，ψ为太阳方向与观测方向之间的相对方位角，即两个方向方位角之差；
得到三个辅助方位角后即可计算单叶片反射率与透射率相关的乘子，用于计算双向散射系数ω；
③计算双向散射系数ω(θl)：
在得到辅助方位角后，配合S2中已计算得到的单叶片反射率和透射率计算得到双向散射系数，计算公式为：
其中，ρ、τ、θl分别表示叶片的反射率、透射率及此时相应的叶倾角；令 S2中的ρ＝Rα(N)；τ＝Tα(N)。
S3.2计算加入叶片反射率和透射率后的散射系数：
①漫辐射E‐和E+的后向散射系数计算公式为：
②漫辐射E-和E+的前向散射系数计算公式为：
③太阳直射辐射ES的后向散射系数计算公式为：
④太阳直射辐射ES的前向散射系数计算公式为：
⑤漫辐射E-和E+的衰减系数计算公式为：
⑥观测方向辐射E0的后向散射系数计算公式为：
⑦观测方向辐射E0的前向散射系数计算公式为：
S3.3计算一般叶倾角分布概率(叶型分布要素)：
计算平均叶倾角分布下的一般叶倾角分布：叶倾角是指水平面法线方向和叶片法线方向夹角；植被冠层的不同类型对应于不同的植被叶片倾斜角度，同时对应于两个叶分布参数LIDFa和LIDFb。
在非椭球形分布的情况下，可根据平均叶倾角计算一般叶倾角分布概率，计算过程通过以下迭代完成：
x＝2θ
y＝LIDFa·sin(x)+0.5LIDFb·sin(2x)
dx＝0.5(y‐x+2θ)
x＝x+dx
直到|dx|＜t；
则F(θ)＝2(y+θ)/π。
其中，θ为平均叶倾角的离散值，F(θ)为累积叶倾角，即一般叶倾角分布概率。
对于椭球形分布，参数LIDFa表示的是分布角度，为30度，参数LIDFb为 0，在求一般叶倾角分布概率时，可利用Campbell(1986)所作的叶倾角密度函数求得；
S3.4确定整个冠层的模型参数；
对于任何冠层来说，叶片倾角不会是固定的，是一个从0到90不断增加的连续过程，因此确定整个冠层的模型参数时，计算方法为将叶倾角为θl时的一般叶倾角分布概率F(θl)与模型参数乘积后再相加，计算公式为：
Z＝∑F(θl)Z(θl)
其中，F(θl)为叶倾角θl对应的一般叶倾角分布概率，Z(θl)为单叶片的参数， Z是经过叶倾角修正模型Z＝∑F(θl)Z(θl)加以修正后的整个冠层的参数；，Z(θl) 指代步骤S3.1和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具体包括以下步骤：
S4.1：计算叶面积参数；叶面积参数是评价植被冠层反射率的重要指标，其参数计算公式为：
τss＝e-klai
τoo＝e-Klai
ρdd＝(emlai-e-mlai)/(h1emlai-h2e-mlai)
τdd＝(h1-h2)/(h1emlai-h2e-mlai)
其中，K和k分别表示观测方向辐射E0的消光系数和太阳直射辐射Es的消光系数；lai表示叶面积指数，为模型输入参数；
S4.2：计算加入热点效应后的参数，计算公式如下：
ρsd＝CS(1-τssτdd)-DSρdd
τsd＝DS(τss-τdd)-CSτssρdd
ρdo＝Co(1-τooτdd)-Doρdd
τdo＝Do(τoo-τdd)-Coτooρdd
ρso＝Ho(1-τssτoo)-Coτsdτoo-Doρsd
其中，各个中间变量的计算方法为：
h1＝(a+m)/σ
h2＝(a-m)/σ＝1/h1
CS＝[s′(k-a)-sσ]/(k2-m2)
Co＝[v(K-a)-uσ]/(K2-m2)
Ds＝[-s(k+a)-s′σ]/(k2-m2)
Do＝[-u(K+a)-vσ]/(K2-m2)
Hs＝(uCS+vDs+w)/(K+k)
Ho＝(sCo+s′Do+w)/(K+k)
其中，τssoo为热点效应修正参数；由于太阳辐射的问题，特别是稀疏植被区域的土壤和植被的温差相当大，这与其物理表面和气象特点相关，因此需要一个温度修正的过程，同时也可能提供额外的植被冠层信息。热点效应的修正过程为：
根据2/(k+K)效应进行热点纠正，首先根据给定的热点值计算参数：
给定alf的初始值为alf＝106，
1)若热点值q大于0，则认为是有热点效应影响，此时 其中，θs为太阳天顶角、θo为观测天顶角、为太阳方位角；(若α计算结果大于200则其alf值为200)；若alf的计算结果为0，则按以下公式计算出无阴影纯热点效应影响下的参数τssoo和sumint：
τssoo＝τss
若alf的计算结果不为0，则认为是无热点效应影响，按2)中的方法计算τssoo和sumint；
2)若热点值q为0，则认为是无热点效应影响，利用循环相加的方法计算出无热点效应影响下的参数τssoo和sumint：具体步骤为:
第一步、参数初始化：
第二步、通过以下迭代计算τ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：引入土壤反射作用：
土壤位于整个冠层反射***的最底部，电磁波穿透冠层后会射到土壤地面，经反射后再次反射回冠层顶部，形成冠层反射率的一部分，计算方案为：
dn＝1－rs*ρdd
S4.5：计算双半球反射率因子：
S4.6：计算太阳直射方向半球反射率：
ρsd2＝ρdd+(τsd+τss)*rs*τdd/dn
S4.7：计算观测方向半球反射率：
ρdo2＝ρdo+(τdo+τoo)*rs*τdd/dn
S4.8：计算双直射反射率：
ρsod2＝ρso+((τss+τsd)*τdo+(τsd+τss*rs*ρdd)*τoo)*rs/dn
ρso3＝ρso2+ρso+τssoo*rs
其中，rs是土壤反射率。
所述步骤S5具体包括以下步骤：
首先，根据数据库资料可以求得太阳直射光线与大气散射光线的基本辐射量，分别记为Es和Ed；
然后，计算表示大气中散射辐射占总辐射比例的系数skyl：
skyl＝0.847-1.61*sin(90-θs)+1.04*sin(90-θs)2
最后，计算观测方向测定到的辐射量与入射的辐射量的比值，得到冠层反射率，公式为：
一种阔叶植被冠层反射率的计算模型，公式如下：
其中，R为冠层反射率，Es和Ed分别为太阳直射光线与大气散射光线的基本辐射量，skyl为表示大气中散射辐射占总辐射比例的系数；ρdo2和ρso3分别为观测方向半球反射率和直射反射率；模型参数根据上述阔叶植被冠层反射率的计算方法求解。
有益效果：
本发明以最新的PROSPECT阔叶模型及SAIL冠层模型研究成果为基础，将模拟阔叶叶片光谱信息的PROSPECT模型与模拟冠层光谱信息的SAIL模型进行耦合，充分利用可获取的参数省去叶片光谱模拟这一过程，取消植被冠层反射率模拟过程中叶片反射率及透射率输入的过程，将叶片模型模拟算法加入冠层 SAIL模型完成耦合过程，直接由植被理化参数模拟冠层光谱。利用物理模型最新研究成果进行耦合可大大提高物理模型模拟光谱的精度，提高模型业务化动作能力，同时有效简化植被冠层光谱信息模拟过程中的参数获取问题并优化算法，加速计算过程，提高植被信息的反演能力。
附图说明
图1为本发明原理图；
图2为本发明输入的模型参数；
图3为本发明实验结果
具体实施方式
以下结合附图对本发明进行进一步说明。
本发明采用热点效应改进后的SAILH冠层模型为基础，首先将给定的参数进行分类，分为叶片参数和冠层参数两部分，之后利用叶片的理化参数，计算叶片的基础吸收作用，再结合无界面吸收作用下的透射率计算导出单层平板的反射率与透射率，最后利用分层理论导出给定结构参数下叶片的光谱信息；同时，对于冠层参数首先利用给定的叶倾角分布模型参数计算叶型分布要素，再结合计算出的几何要素计算冠层模型的消光系数及散射系数并与叶型分布模型结合，经过循环过程得到完整叶片分布下的成套冠层参数，最后，利用求得的两个辐射参数及天空散射系数求定冠层的反射率。
本实施案例以高速公路边某颗生长情况良好且冠层光谱容易测定的树为例对本发明的技术方案进行进一步说明。如附图所示，本发明所涉及的改进型物理模型PROSAILH的计算方法如下：
经过测定，该树的参数如下表：
叶绿素含量 45.80 叶倾角分布类型 倾斜型
等效水厚度 0.02 叶面积指数 4.30
干物质含量 0.02 太阳天顶角 30.00
叶片结构参数 1.30 观测天顶角 10.00
热点 8.00 相对方位角 0.00
选择实验波长为400nm-2399nm
S1参数分类：在得到植被的相关参数后分为两部分，其中叶片参数包括叶片的色素含量、等效水厚度、干物质含量及最大入射角(一般为59度)；冠层参数包括叶倾角分布参数、叶面积指数、太阳天顶角、观测天顶角、相对方位角及热点和土壤的背景反射率。
S2：叶片光谱信息解算：在取得相应参数后，首先需要计算单个叶片的光谱反射率及透射率，计算方案如下：
(1)吸收系数K计算：吸收系数是叶片各种化学物质含量的线性表达，计算公式为：
(2)无界面吸收作用下透射率计算：自然光线可看作是非极化光，反射后发生完全极化，透射后折射光线发生部分极化。根据Snell-Descarts定律可知道两种介质接触面上入射角与折射角及两个折射率之间的关系，再根据 Schanda给出的电磁波穿过两种介质的透射率即可计算出光线以α立体角由折射率为1的空气入射折射率为n的叶片的平均透射率tav(α,1,n),此时赋予其代码t12,根据Stern的研究成果，光线从折射率为n的叶片射入折射率为1的空气透射率t21与t12之间存在关系式t21＝n－2t12，据此计算出界面处的透射率。
(3)光线在各向同性介质中的传输：对于自然光，其透过平板介质的总透射率τ为整个半球2π空间内的积分，与吸收系数K及叶片厚度D相关，计算公式为：
(4)单层平板光谱模拟：单层平板的光谱是整个叶片光谱模拟的核心，主要包括单层平板的反射率及透射率，其计算公式为：
(5)单片叶反射率及透射率模拟：Stokes对光线穿透有限层平板后的光学现象进行了深入研究，得出了光线穿过N层相同反射率和透射率的平板后的总体反射率和透射率，此时应用其理论可得出叶片光谱信息的模拟公式：
其中：
S3：计算消光系数及散射系数
(6)几何要素计算：
(7)叶型分布要素计算：
平均叶倾角分布下的一般叶倾角分布：植被冠层的不同类型对应于不同的植被叶片倾斜角度，同时对应于两个叶分布参数LIDFa和LIDFb。在非椭球形分布的情况下，可根据平均叶倾角计算一般倾角，计算过程可用一个简单的迭代完成，伪代码为：
x＝2θ
y＝LIDFa·sin(x)+0.5LIDFb·sin(2x)
dx＝0.5(y‐x+2θ)
x＝x+dx
直到|dx|＜t；
则F(θ)＝2(y+θ)/π
其中θ为平均叶倾角的离散值，F(θ)为累积叶倾角。
对于椭球形分布，参数LIDFa表示的是分布角度，为30度，参数LIDFb 为0，在求一般角度分布时，可利用Campbell(1986)所作的叶倾角密度函数求得。
(8)消光及散射系数计算：
阴影补偿：
①与消光和散射相关的几何因子计算：首先给叶倾角分布加权并离散化，得到一组离散化的中心角度值叶倾角分布。之后，每个叶倾角对应计算一组消光与散射因子，单次计算方案如下：
首先求临界角βs和βo，计算公式为：
当确定计算公式中分母不为0且cosβ的计算结果小于1时直接计算出两个临界角的值；当cosβ的计算结果大于1时，两个临介角均等于π；其中cosβ泛指βs和βo的余弦；
计算出两个临时界角后就可以计算太阳直射方向与观测方向的消光系数，计算公式分别为：
其中，L′＝lai/h，lai为叶面积指数，h为冠层高度；
②用于计算双向散射系数W的辅助方位角β1、β2、β3计算：辅助方位角的计算主要取决于两个临界角的相对不等关系，取值方法如下：
If： β1 β2 β3
ψ≤|β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)。
模型参数解算
(9)散射系数的计算(叶片反射率和透射率的加入)：
①漫辐射E‐和E+的后向散射系数计算公式为：
②漫辐射E-和E+的前向散射系数计算公式为：
③太阳直射辐射ES的后向散射系数计算公式为：
④太阳直射辐射ES的前向散射系数计算公式为：
⑤漫辐射E-和E+的衰减系数计算公式为：
⑥观测方向E0的后向散射系数计算公式为：
⑦观测方向E0的前向散射系数计算公式为：
对于任何冠层来说，叶片倾角不会是固定的，是一个从0到90不断增加的连续过程，因此确定整个冠层的模型参数时，计算方法为叶倾角分布函数作为概率函数F(θ)与模型参数乘积后再相加，计算公式为：
Z＝∑F(θl)Z(θl)
其中，F(θl)为叶倾角θl对应的一般叶倾角分布概率，Z(θl)为单叶片的参数， Z是经过叶倾角修正模型Z＝∑F(θl)Z(θl)加以修正后的整个冠层的参数；，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＝(emlai-e-mlai)/(h1emlai-h2e-mlai)
τdd＝(h1-h2)/(h1emlai-h2e-mlai)
(2)解的参数计算：加入热点效应后的参数较复杂，计算公式如下：
ρsd＝CS(1-τssτdd)-DSρdd
τsd＝DS(τss-τdd)-CSτssρdd
ρdo＝Co(1-τooτdd)-Doρdd
τdo＝Do(τoo-τdd)-Coτooρdd
ρso＝Ho(1-τssτoo)-Coτsdτoo-Doρsd
其中各个变量的计算方法为：
h1＝(a+m)/σ
h2＝(a-m)/σ＝1/h1
CS＝[s′(k-a)-sσ]/(k2-m2)
Co＝[v(K-a)-uσ]/(K2-m2)
Ds＝[-s(k+a)-s′σ]/(k2-m2)
Do＝[-u(K+a)-vσ]/(K2-m2)
Hs＝(uCS+vDs+w)/(K+k)
Ho＝(sCo+s′Do+w)/(K+k)
(3)热点效应的修正：由于太阳辐射的问题，特别是稀疏植被区域的土壤和植被的温差相当大，这与其物理表面和气象特点相关，因此需要一个温度修正的过程，同时也可能提供额外的植被冠层信息。热点效应的修正过程为：
根据2/(k+K)效应进行热点纠正，首先根据给定的热点值计算参数：
给定alf的初始值为alf＝106，
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－rs*ρdd
(6)双半球反射率因子计算：
(7)太阳直射方向半球反射率计算：
ρsd2＝ρdd+(τsd+τss)*rs*τdd/(1－rs*ρdd)
(8)观测方向半球反射率计算：
ρdo2＝ρdo+(τdo+τoo)*rs*τdd/(1－rs*ρdd)
(9)双直射反射率计算：
ρsod2＝ρso+((τss+τsd)*τdo+(τsd+τss*rs*ρdd)*τoo)*rs/(1－rs*ρdd)
ρso3＝ρso2+ρso+τssoo*rs
S5：计算冠层反射率：
根据数据库资料可以将太阳直射光线与大气散射光线的基本辐射量求得，分别为ES和Ed，其次计算大气散射系数skyl参数：
sskyl＝0.847-1.61*sin(90-θs)+1.04*sin(90-θs)2
冠层反射率的计算方法为观测方向测定到的辐射量与入射的辐射量的比值，公式为：
计算出的植被反射率如图3所。
本次选取的高速公路路域植被理化参数及对比所用光谱均为实测值，模型模拟结果利用RMSE确定可用性及精度问题。
经检测，该模型本例模拟的光谱信息均方根误差(RMSE)为0.15706，模拟算法模拟结果可用且精度较高。
所描述的实例是本发明的一部分实例，而不是全部实例，不能理解为对本发明的限制。基于本发明中的实例，本领域普通技术人员在没有做出创新性劳动的前提下所获得的所有其他实施方式都属于本发明的保护范围。

Claims (3)

1.一种阔叶植被冠层反射率的计算方法，其特征在于，包括以下步骤：
S1：参数识别；
输入模型参数；并将输入的参数初步分为三大类：叶片参数、土壤参数和冠层参数；
S2：根据步骤S1中的叶片参数输入PROSPECT模型进行单片叶光谱模拟，计算单个叶片的光谱反射率及透射率；
S3：根据步骤S1中的土壤参数和冠层参数、步骤S2得到的单个叶片的光谱反射率及透射率计算消光系数及散射系数；
S4：将步骤S3计算得到的消光系数及散射系数输入SAIL模型，计算冠层的相关反射因子和反射率；
S5：计算冠层反射率；
所述步骤S2具体包括以下步骤：
S2.1：计算吸收系数K：
<mrow> <mi>K</mi> <mo>=</mo> <mo>&amp;Sigma;</mo> <mfrac> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>*</mo> <mi>k</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mrow> <mi>N</mi> </mfrac> </mrow>
其中，c(i)为叶片中组分i的浓度；k(i)为组分i的特定吸收系数；N为叶片结构参数，表示叶片分层数；N采用以下经验方程计算得到：
N＝(0.9*SLA+0.025)/(SLA-0.1)
其中，SLA为比叶面积，指每单位干重的叶面积；
S2.2：计算界面处的透射率：
计算出光线以立体角α由折射率为1的空气入射折射率为n的叶片的平均透射率tav(α,1,n)，令t12＝tav(α,1,n)；
S2.3：计算光线透过平板介质的总透射率τ1；
<mrow> <msub> <mi>&amp;tau;</mi> <mn>1</mn> </msub> <mo>=</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>K</mi> <mi>D</mi> <mo>)</mo> </mrow> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mi>K</mi> <mi>D</mi> <mo>)</mo> </mrow> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>K</mi> <mi>D</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msubsup> <mo>&amp;Integral;</mo> <mrow> <mi>K</mi> <mi>D</mi> </mrow> <mi>&amp;infin;</mi> </msubsup> <msup> <mi>t</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>exp</mi> <mrow> <mo>-</mo> <mi>t</mi> </mrow> </msup> <mi>d</mi> <mi>t</mi> </mrow>
其中，t为中间变量，K为吸收系数，D为叶片厚度；
S2.4：计算单层平板的反射率Rα(1)及透射率Tα(1)：
<mrow> <msub> <mi>R</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>t</mi> <mn>12</mn> </msub> <mo>+</mo> <mfrac> <mrow> <msup> <msub> <mi>&amp;tau;</mi> <mn>1</mn> </msub> <mn>2</mn> </msup> <msubsup> <mi>t</mi> <mn>12</mn> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>-</mo> <msub> <mi>t</mi> <mn>12</mn> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mi>n</mi> <mn>4</mn> </msup> <mo>-</mo> <msup> <msub> <mi>&amp;tau;</mi> <mn>1</mn> </msub> <mn>2</mn> </msup> <msup> <mrow> <mo>(</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>-</mo> <msub> <mi>t</mi> <mn>12</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mfrac> </mrow>
<mrow> <msub> <mi>T</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msup> <mi>n</mi> <mn>2</mn> </msup> <msup> <msub> <mi>&amp;tau;</mi> <mn>1</mn> </msub> <mn>2</mn> </msup> <msubsup> <mi>t</mi> <mn>12</mn> <mn>2</mn> </msubsup> </mrow> <mrow> <msup> <mi>n</mi> <mn>4</mn> </msup> <mo>-</mo> <msup> <msub> <mi>&amp;tau;</mi> <mn>1</mn> </msub> <mn>2</mn> </msup> <msup> <mrow> <mo>(</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>-</mo> <msub> <mi>t</mi> <mn>12</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mfrac> </mrow>
其中，n为叶片的折射率；
S2.5：计算单片叶反射率Rα(N)及透射率Tα(N)，即光线穿过N层相同反射率和透射率的平板后的总体反射率和透射率：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>R</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mi>N</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>R</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msup> <mi>ab</mi> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <msup> <mi>a</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>b</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>N</mi> </mrow> </msup> <mo>)</mo> </mrow> <mo>+</mo> <mrow> <mo>(</mo> <msub> <mi>T</mi> <mi>&amp;alpha;</mi> </msub> <mo>(</mo> <mn>1</mn> <mo>)</mo> <msub> <mi>T</mi> <mn>90</mn> </msub> <mo>(</mo> <mn>1</mn> <mo>)</mo> <mo>-</mo> <msub> <mi>R</mi> <mi>&amp;alpha;</mi> </msub> <mo>(</mo> <mn>1</mn> <mo>)</mo> <msub> <mi>R</mi> <mn>90</mn> </msub> <mo>(</mo> <mn>1</mn> <mo>)</mo> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msup> <mi>b</mi> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <msup> <mi>b</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>N</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mi>ab</mi> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <msup> <mi>a</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>b</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>N</mi> </mrow> </msup> <mo>-</mo> <msub> <mi>R</mi> <mn>90</mn> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msup> <mi>b</mi> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <msup> <mi>b</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>N</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mrow> <msub> <mi>T</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mi>N</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>T</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>a</mi> <mo>-</mo> <msup> <mi>a</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mi>ab</mi> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <msup> <mi>a</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>b</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>N</mi> </mrow> </msup> <mo>-</mo> <msub> <mi>R</mi> <mn>90</mn> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msup> <mi>b</mi> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <msup> <mi>b</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>N</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
其中：
<mrow> <mi>&amp;delta;</mi> <mo>=</mo> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>R</mi> <mn>90</mn> </msub> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <msub> <mi>T</mi> <mn>90</mn> </msub> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mn>4</mn> <msub> <mi>T</mi> <mn>90</mn> </msub> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mo>;</mo> </mrow>
所述步骤S3具体包括以下步骤：
S3.1阴影补偿：
①计算与消光和散射相关的几何因子：
首先根据一般叶倾角分布概率加权并离散化，得到一组离散化的叶倾角(即水平面法线方向和叶片法线方向夹角)；之后，计算每个叶倾角对应的消光与散射因子；对于叶倾角θl，计算方法如下：
首先，求临界角βs和βo，计算公式为：
<mrow> <msub> <mi>&amp;beta;</mi> <mi>s</mi> </msub> <mo>=</mo> <mi>arccos</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>cos&amp;theta;</mi> <mi>l</mi> </msub> <msub> <mi>cos&amp;theta;</mi> <mi>s</mi> </msub> </mrow> <mrow> <msub> <mi>sin&amp;theta;</mi> <mi>l</mi> </msub> <msub> <mi>sin&amp;theta;</mi> <mi>s</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> </mrow>
<mrow> <msub> <mi>&amp;beta;</mi> <mi>o</mi> </msub> <mo>=</mo> <mi>arccos</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>cos&amp;theta;</mi> <mi>l</mi> </msub> <msub> <mi>cos&amp;theta;</mi> <mi>o</mi> </msub> </mrow> <mrow> <msub> <mi>sin&amp;theta;</mi> <mi>l</mi> </msub> <msub> <mi>sin&amp;theta;</mi> <mi>o</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> </mrow>
其中，θs为太阳天顶角、θo为观测天顶角；
当确定计算公式中分母不为0且cosβ的计算结果小于1时直接计算出两个临界角的值；当cosβ的计算结果等于1时，两个临介角均等于π；其中cosβ泛指βs和βo的余弦；
然后，计算叶倾角为θl的单个叶片太阳直射方向与观测方向的消光系数，计算公式分别为：
<mrow> <mi>k</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>2</mn> <mi>&amp;pi;</mi> </mfrac> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <mo>&amp;lsqb;</mo> <mrow> <mo>(</mo> <msub> <mi>&amp;beta;</mi> <mi>s</mi> </msub> <mo>-</mo> <mfrac> <mi>&amp;pi;</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <msub> <mi>cos&amp;theta;</mi> <mi>l</mi> </msub> <mo>+</mo> <msub> <mi>sin&amp;beta;</mi> <mi>s</mi> </msub> <msub> <mi>tan&amp;theta;</mi> <mi>s</mi> </msub> <msub> <mi>sin&amp;theta;</mi> <mi>l</mi> </msub> <mo>&amp;rsqb;</mo> </mrow>
<mrow> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>2</mn> <mi>&amp;pi;</mi> </mfrac> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <mo>&amp;lsqb;</mo> <mrow> <mo>(</mo> <msub> <mi>&amp;beta;</mi> <mi>o</mi> </msub> <mo>-</mo> <mfrac> <mi>&amp;pi;</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <msub> <mi>cos&amp;theta;</mi> <mi>l</mi> </msub> <mo>+</mo> <msub> <mi>sin&amp;beta;</mi> <mi>o</mi> </msub> <msub> <mi>tan&amp;theta;</mi> <mi>o</mi> </msub> <msub> <mi>sin&amp;theta;</mi> <mi>l</mi> </msub> <mo>&amp;rsqb;</mo> </mrow>
其中，L′＝lai/h，lai为叶面积指数，h为冠层高度；
②计算辅助方位角β1、β2、β3：取值方法如下：
If： β1β β2β β3β ψ≤|β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 ψ
其中，ψψ为太阳方向与观测方向之间的相对方位角，即两个方向方位角之差；
③计算双向散射系数ω(θl)：
在得到辅助方位角后，配合S2中已计算得到的单叶片反射率Rα(N)及透射率Tα(N)和透射率计算得到双向散射系数，计算公式为：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>&amp;omega;</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mn>2</mn> <mi>&amp;pi;</mi> </mrow> </mfrac> <mo>{</mo> <mo>&amp;lsqb;</mo> <mi>&amp;pi;</mi> <mi>&amp;rho;</mi> <mo>-</mo> <msub> <mi>&amp;beta;</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>+</mo> <mi>p</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mrow> <mo>(</mo> <mn>2</mn> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>+</mo> <msup> <mi>sin</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <msub> <mi>tan&amp;theta;</mi> <mi>s</mi> </msub> <msub> <mi>tan&amp;theta;</mi> <mi>o</mi> </msub> <mi>cos</mi> <mi>&amp;Psi;</mi> <mo>)</mo> </mrow> <mo>+</mo> <mo>(</mo> <mi>&amp;rho;</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mi>&amp;tau;</mi> <mo>)</mo> <msub> <mi>sin&amp;beta;</mi> <mn>2</mn> </msub> <mo>&amp;lsqb;</mo> <mfrac> <mrow> <mn>2</mn> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> </mrow> <mrow> <msub> <mi>cos&amp;beta;</mi> <mi>s</mi> </msub> <msub> <mi>cos&amp;beta;</mi> <mi>o</mi> </msub> </mrow> </mfrac> <mo>+</mo> <msub> <mi>cos&amp;beta;</mi> <mn>1</mn> </msub> <msub> <mi>cos&amp;beta;</mi> <mn>3</mn> </msub> <msup> <mi>sin</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <msub> <mi>tan&amp;theta;</mi> <mi>s</mi> </msub> <msub> <mi>tan&amp;theta;</mi> <mi>o</mi> </msub> <mo>&amp;rsqb;</mo> <mo>}</mo> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，ρ、τ、θl分别表示叶片的反射率、透射率及此时相应的叶倾角；ρ＝Rα(N)；τ＝Tα(N)；
S3.2计算加入叶片反射率和透射率后的散射系数：
①漫辐射E-和E+的后向散射系数计算公式为：
<mrow> <mi>&amp;sigma;</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&amp;tau;</mi> <mo>+</mo> <mi>&amp;rho;</mi> </mrow> <mn>2</mn> </mfrac> <mo>+</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>-</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> </mrow>
②漫辐射E-和E+的前向散射系数计算公式为：
<mrow> <msup> <mi>&amp;sigma;</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&amp;tau;</mi> <mo>+</mo> <mi>&amp;rho;</mi> </mrow> <mn>2</mn> </mfrac> <mo>-</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>-</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> </mrow>
③太阳直射辐射ES的后向散射系数计算公式为：
<mrow> <mi>s</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>+</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <mi>k</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>-</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> </mrow>
④太阳直射辐射ES的前向散射系数计算公式为：
<mrow> <msup> <mi>s</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>+</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <mi>k</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>-</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> </mrow>
⑤漫辐射E-和E+的衰减系数计算公式为：
<mrow> <mi>a</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>+</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <mo>+</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>-</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> </mrow>
⑥观测方向辐射E0的后向散射系数计算公式为：
<mrow> <mi>v</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>+</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>-</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> </mrow>
⑦观测方向辐射E0的前向散射系数计算公式为：
<mrow> <mi>u</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>+</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mo>-</mo> <mi>&amp;tau;</mi> </mrow> <mn>2</mn> </mfrac> <msup> <mi>L</mi> <mo>&amp;prime;</mo> </msup> <msup> <mi>cos</mi> <mn>2</mn> </msup> <msub> <mi>&amp;theta;</mi> <mi>l</mi> </msub> </mrow>
S3.3首先针对①中离散化得到的每一个叶倾角θ，计算其对应的一般叶倾角分布概率F(θ)：
在非椭球形分布的情况下，根据平均叶倾角计算一般叶倾角分布概率，计算过程通过以下迭代完成：
x＝2θ
y＝LIDFa·sin(x)+0.5LIDFb·sin(2x)
dx＝0.5(y-x+2θ)
x＝x+dx
直到|dx|<t；
则F(θ)＝2(y+θ)/π；
其中，LIDFa和LIDFb为叶分布参数；θ为平均叶倾角的离散值，F(θ)为累积叶倾角，即一般叶倾角分布概率；
对于椭球形分布，利用Campbell叶倾角密度函数求一般叶倾角分布概率F(θ)；
S3.4确定整个冠层的模型参数；
将①中离散化得到的各个叶倾角对应的叶倾角分布概率分别与模型参数相乘后再相加，得到新的模型参数，计算公式为：
Z＝∑F(θl)Z(θl)
其中，F(θl)为叶倾角θl对应的一般叶倾角分布概率，Z(θl)指代步骤S3.1和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中的任意一个。
2.根据权利要求1所述的阔叶植被冠层反射率的计算方法，其特征在于，所述步骤S4具体包括以下步骤：
S4.1：计算叶面积参数：
τss＝e-klai
τoo＝e-Klai
ρdd＝(emlai-e-mlai)/(h1emlai-h2e-mlai)
τdd＝(h1-h2)/(h1emlai-h2e-mlai)
其中，lai表示叶面积指数，为模型输入参数；
S4.2：计算加入热点效应后的参数，计算公式如下：
ρsd＝CS(1-τssτdd)-DSρdd
τsd＝DS(τss-τdd)-CSτssρdd
ρdo＝Co(1-τooτdd)-Doρdd
τdo＝Do(τoo-τdd)-Coτooρdd
ρso＝Ho(1-τssτoo)-Coτsdτoo-Doρsd
其中，各个中间变量的计算方法为：
<mrow> <mi>m</mi> <mo>=</mo> <msqrt> <mrow> <msup> <mi>a</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>&amp;sigma;</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mrow>
h1＝(a+m)/σ
h2＝(a-m)/σ＝1/h1
CS＝[s′(k-a)-sσ]/(k2-m2)
Co＝[v(K-a)-uσ]/(K2-m2)
Ds＝[-s(k+a)-s′σ]/(k2-m2)
Do＝[-u(K+a)-vσ]/(K2-m2)
Hs＝(uCS+vDs+w)/(K+k)
Ho＝(sCo+s′Do+w)/(K+k)
其中，τssoo和sumint为热点效应修正参数；
S4.3：计算双向反射率：
双向反射率包括两部分，一部分带热点效应的单次散射影响，另一部分未带热点效应的多次散射影响，两部分求和为最终对冠层顶部的反射率影响的计算结果：
ρso2＝ρso+w*lai*sumint
S4.4：引入土壤反射作用：
dn＝1－rs*ρdd
S4.5：计算双半球反射率因子：
<mrow> <msub> <mi>&amp;rho;</mi> <mrow> <mi>dd</mi> <mn>2</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;rho;</mi> <mi>dd</mi> </msub> <mo>+</mo> <mrow> <mo>(</mo> <msubsup> <mi>&amp;tau;</mi> <mi>dd</mi> <mn>2</mn> </msubsup> <mo>*</mo> <msub> <mi>r</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>d</mi> <mi>n</mi> </msub> </mrow>
S4.6：计算太阳直射方向半球反射率：
ρsd2＝ρdd+(τsd+τss)*rs*τdd/dn
S4.7：计算观测方向半球反射率：
ρdo2＝ρdo+(τdo+τoo)*rs*τdd/dn
S4.8：计算双直射反射率：
ρsod2＝ρso+((τss+τsd)*τdo+(τsd+τss*rs*ρdd)*τoo)*rs/dn
ρso3＝ρso2+ρso+τssoo*rs
其中，rs是土壤反射率。
3.根据权利要求2所述的阔叶植被冠层反射率的计算方法，其特征在于，所述步骤S5具体包括以下步骤：
首先，根据数据库资料求得太阳直射光线与大气散射光线的基本辐射量，分别记为Es和Ed；
然后，计算表示大气中散射辐射占总辐射比例的系数skyl：
skyl＝0.847-1.61*sin(90-θs)+1.04*sin(90-θs)2
最后，计算观测方向测定到的辐射量与入射的辐射量的比值，得到冠层反射率，公式为：
<mrow> <mi>R</mi> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&amp;rho;</mi> <mrow> <mi>d</mi> <mi>o</mi> <mn>2</mn> </mrow> </msub> <mo>*</mo> <mi>s</mi> <mi>k</mi> <mi>y</mi> <mi>l</mi> <mo>*</mo> <msub> <mi>E</mi> <mi>d</mi> </msub> <mo>+</mo> <msub> <mi>&amp;rho;</mi> <mrow> <mi>s</mi> <mi>o</mi> <mn>3</mn> </mrow> </msub> <mo>*</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>s</mi> <mi>k</mi> <mi>y</mi> <mi>l</mi> <mo>)</mo> </mrow> <mo>*</mo> <msub> <mi>E</mi> <mi>s</mi> </msub> </mrow> <mrow> <mi>s</mi> <mi>k</mi> <mi>y</mi> <mi>l</mi> <mo>*</mo> <msub> <mi>E</mi> <mi>d</mi> </msub> <mo>+</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>s</mi> <mi>k</mi> <mi>y</mi> <mn>1</mn> <mo>)</mo> </mrow> <mo>*</mo> <msub> <mi>E</mi> <mi>s</mi> </msub> </mrow> </mfrac> <mo>.</mo> </mrow>
CN201710141821.0A 2017-03-10 2017-03-10 一种阔叶植被冠层反射率的计算方法 Active CN106909750B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710141821.0A CN106909750B (zh) 2017-03-10 2017-03-10 一种阔叶植被冠层反射率的计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710141821.0A CN106909750B (zh) 2017-03-10 2017-03-10 一种阔叶植被冠层反射率的计算方法

Publications (2)

Publication Number Publication Date
CN106909750A CN106909750A (zh) 2017-06-30
CN106909750B true CN106909750B (zh) 2018-02-27

Family

ID=59186335

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710141821.0A Active CN106909750B (zh) 2017-03-10 2017-03-10 一种阔叶植被冠层反射率的计算方法

Country Status (1)

Country Link
CN (1) CN106909750B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109033562B (zh) * 2018-07-05 2020-08-11 浙江大学 一种卷曲状态下叶片二向反射值的计算方法
CN111007024B (zh) * 2019-12-25 2021-01-26 武汉大学 一种适用于氧气a带的云反射比快速确定方法
CN112254820B (zh) * 2020-10-14 2021-04-27 中国科学院空天信息创新研究院 一种离散森林场景热红外辐射传输模拟方法
CN112666121B (zh) * 2020-12-17 2024-04-05 淮阴师范学院 基于红外光谱的植被与非植被识别方法
CN114528728B (zh) * 2022-01-20 2024-05-28 北京航空航天大学 一种行播水生植被冠层反射的蒙特卡罗计算机模拟方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103630495B (zh) * 2013-11-13 2015-08-12 北京航空航天大学 一种水生植被-大气耦合辐射传输模型
CN103632040B (zh) * 2013-11-14 2016-05-04 北京航空航天大学 一种通用的水生植被辐射传输模型

Also Published As

Publication number Publication date
CN106909750A (zh) 2017-06-30

Similar Documents

Publication Publication Date Title
CN106909750B (zh) 一种阔叶植被冠层反射率的计算方法
Li et al. Multi-LUTs method for canopy nitrogen density estimation in winter wheat by field and UAV hyperspectral
CN106874621B (zh) 一种针叶植被冠层反射率计算方法
CN107437267B (zh) 植被区高光谱图像模拟方法
CN108009392B (zh) 一种浓密植被地表的遥感反射率模型构建及标定应用方法
CN104483271A (zh) 光学反射模型与微波散射模型协同的森林生物量反演方法
CN114460013B (zh) 滨海湿地植被地上生物量gan模型自学习遥感反演方法
CN114581791A (zh) 基于modis数据大气水汽含量反演方法及***
CN108760662A (zh) 一种卫星遥感大气臭氧廓线反演算法
Lee et al. Development of land surface albedo algorithm for the GK-2A/AMI instrument
CN106202971A (zh) 基于folium模型叶片色素遥感反演方法
Lin et al. Developing a two-step algorithm to estimate the leaf area index of forests with complex structures based on CHRIS/PROBA data
CN111220552B (zh) 考虑光照方向叶片辐射传输模型的叶绿素高光谱反演方法
CN112924401B (zh) 一种植被冠层叶绿素含量半经验反演方法
Zhao et al. Variations and driving mechanisms of desertification in the southeast section of the China-Mongolia-Russia Economic Zone
Zeng et al. Dynamic monitoring of plant cover and soil erosion using remote sensing, mathematical modeling, computer simulation and GIS techniques
Bouin et al. Using scintillometry to estimate sensible heat fluxes over water: First insights
Nur et al. Comparison of soil moisture content retrieval models utilizing hyperspectral goniometer data and hyperspectral imagery from an unmanned aerial system
CN116223452A (zh) 基于卷曲植物叶片辐射传输模型的叶绿素高光谱反演方法
Wang et al. Study on Remote Sensing of Water Depths Based on BP Artificial Neural Network.
Edalati et al. Estimating and modeling monthly mean daily global solar radiation on horizontal surfaces using artificial neural networks
CN105259145A (zh) 一种同时遥感岛礁水下地形和地物的方法
Mridha et al. Comparative evaluation of inversion approaches of the radiative transfer model for estimation of crop biophysical parameters
CN107274460B (zh) 一种全谱段高光谱图像模拟方法及装置
Guo et al. New two-step species-level AGB estimation model applied to urban parks

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information

Inventor after: Guo Yunkai

Inventor after: An Guanxing

Inventor after: Liu Haiyang

Inventor after: Jiang Ming

Inventor after: Xie Qiong

Inventor after: Zhou Fengsong

Inventor before: Guo Yunkai

Inventor before: An Guanxing

Inventor before: Xie Qiong

Inventor before: Zhou Fengsong

Inventor before: Li Jian

Inventor before: Zhu Luhong

CB03 Change of inventor or designer information
GR01 Patent grant
GR01 Patent grant