马红章，孙淑怡，刘素美，艾璐，孙根云，孙林

1.中国石油大学 理学院,青岛 266580;

2.中国石油大学 海洋与空间信息学院,青岛 266580;

3.山东科技大学 测绘学院,青岛 266590

1 引 言

植被在地球系统中扮演着重要的角色，是重要的可再生资源（Zhu 等，2019；Zhou 等，2017；马培培等，2019）。植被影响地气系统的能量平衡，是气候、水文和生化循环中最活跃和最具有价值的载体，也是环境评估的敏感指标，因此，植被一直是遥感领域的关注焦点（彭志兴等，2016；Gommes 等，2016；赵燕红 等，2021）。植被冠层空间分布的非均一性导致冠层的二向反射率因子BRF（Bidirectional Reflectance Function）特征变得十分复杂，尤其是随着光学遥感步入亚米时代，基于均一性冠层假设的辐射传输模型已经不再适用，因此，构建异质性冠层三维辐射传输模型对高空间分辨率光学遥感数据的理解和应用具有一定的现实指导意义（徐希孺等，2017；Yan 等，2012；柳钦火 等，2019；Wu 等，2017）。

依托计算机强大的计算能力，基于蒙特卡洛模拟MC（Monte Carlo Method）算法和真实场景模拟的冠层三维辐射传输计算机模型得到了长足发展和普遍应用，例如RAYTRAN 模型（Govaerts等，1998）、DART 模型（Gastellu−Etchegorry 等，2004）、DIANA模型（España 等，1999）以及RAPID模型（Qin 和Gerstl，2000；Huang 等，2013）等。依据算法原理的不同，计算机模型也存在各自的优势与不足（Jiao等，2016）。MC 类模型把光线在冠层内的传输视为众多光子与冠层组分相互作用的统计结果，可实现从光子产生到消亡过程的独立追踪，对光能量在冠层内的传输路径有明确的刻画，因此该类模型对研究冠层对光的散射机理比较有利，但只有样本数足够大时，模型计算精度才能得到保证（李骁尧和黄华国，2020）；基于真实场景的模型对冠层场景刻画更加细致，对冠层叶片组分的空间非均一性分布刻画更逼近现实，但由于模型要求标记每片叶子的空间位置及其取向等参数，导致模型的场景构建非常繁杂，模拟精度也直接受场景刻画逼真度的影响（马红章等，2016）。

本文在融合MC 模型和真实场景模型优势的基础上，利用实际植株统计特征的概率株来构建三维植被冠层，结合蒙特卡洛光线追踪技术和冠层孔隙率的计算，建立了三维冠层辐射传输模型，通过采用‘能量光粒子’代替传统MC 模型中的单光子以及逐步统计粒子能量对观测贡献两种策略来提高模型的稳定性。

2 原理与方法

2.1 冠层构建方法

概率株是大量同类植株几何特征及结构参数的抽象与平均，是刻画冠层的基本单元。概率株被简化为具有中心对称特征的株冠和竖直圆柱形株秆组成的几何体，用株高、最大冠幅半径、冠下高、秆半径和单株叶片总面积等参数来刻画。设定最大冠幅半径为单位1，将株冠等分为10 个水平层，依据植株外形先验知识分别设定各分层冠幅半径的相对比例来刻画株冠的几何边界。株冠内叶片体密度分布的刻画方法是设定株冠内叶片体密度最大值为单位1，依据实际植株叶片垂直向分布的先验知识设定各水平层最大叶片体密度相对比例；将水平层沿径向进行10 等分，依据先验知识设定各等分点处的叶片体密度相对比例，并假定各水平层的叶片体密度径向分布规律相同且不随方位发生改变。通过插值的方法可得到株冠空间内任意点的归一化叶片体密度值（马红章等，2016）。概率株叶倾角分布特征在整个株冠空间内相同。图1展示了叶片体密度在垂直向以及径向均为线性变化的锥形株冠概率株叶片体密度的空间分布。三维冠层是由众多给定位置坐标的概率株组成，若冠层空间内某点A 处于多概率株立体空间重叠位置，则A 点处叶片体密度为各概率株在该点体密度值的累加。

图1 概率株及叶片体密度分布Fig.1 Schematic diagram of probability plant and the distribution of leaf density

2.2 模型构建方法

本文模型将光线视为具有频率和能量双重属性的光粒子，频率属性使模型可同时运算不同频率属性下光粒子在冠层内的传输行为，能量属性可确保光粒子能多次与叶片发生作用直至能量消失或者光粒子逃逸出冠层上边界。概率株位置、叶片角度、光粒子运动方向以及观测几何等参数的定义均基于统一的笛卡尔坐标系FD，一般定义坐标系FD的Z 轴指向天顶，方位角以X 轴为0°起点，逆时针方位增加；若组成冠层的植株分布具有典型的垄行结构，则坐标系FD的X 轴取垂直垄向。光粒子运动垂直向上为0°，垂直向下为180°；叶片的法线方向与Z轴夹角为叶片倾角，在叶片法线为z'轴的叶片坐标系fL下计算光粒子与叶片碰撞前后的运动方向。

2.2.1 叶片散射特性描述

在本文模型中，定义叶片反射率R0和透过率T0分别为光粒子垂直叶片入射时的后向和前向的半球散射比，可通过实验测量或由PROSPECT 模型计算。

叶片对光粒子能量的吸收比则为1 −T0−R0且假定不随入射角度发生变化，前向的散射为朗伯散射即光粒子与叶片碰撞后在前向半球空间各方向出现的概率均等。T随入射角度θ的变化如式（1）：

后向散射为朗伯散射与镜向反射的叠加，两者发生的概率之比取决于叶片的二向散射特性。R随入射角度θ的变化如式（2）：

2.2.2 光粒子运动轨迹计算

模型将冠层水平分层，层厚度足够小满足在层内仅发生单次散射。光粒子在冠层内的传输过程如图2所示。光粒子从冠层上边界发射进入冠层后，由冠层体密度决定光粒子与叶片发生碰撞的概率，用产生机器随机数的方法决定光粒子是否与冠层叶片发生碰撞。若无碰撞，则光粒子沿原方向行进到下一层；若碰撞，由叶片参数R和T决定发生前向散射和后向散射的概率，产生机器随机数决定二者之其一，随之确定光粒子运动方向和行进位置。依次周而复始的迭代计算，直至光粒子能量小于设定阈值（10−5）或者逃逸出冠层上边界。如果光粒子运动到土壤层，则被土壤层吸收部分能量后再朗伯散射返回冠层。如果光粒子从冠层侧面逃逸出冠层，则把光粒子平移到样地坐标原点，继续传输。

图2 光粒子散射过程设计Fig.2 Flow chart of calculation of light particle scattering process

设冠层水平层厚度为d，光粒子由某层的A 点（x（n−1），y（n−1），z（n−1））运动到临近层B点（x（n），y（n），z（n）），与叶片发生作用后移动到C 点（x（n+1），y（n+1），z（n+1）），则C 点坐标的计算分以下3种情况讨论。

（1）无碰撞。若光粒子与叶片没有发生碰撞，则光粒子直接沿原方向行进到C点。

式中，k为方向因子，若碰撞后光粒子运动天顶角<90°，则k=1，否则k=−1。

（2）镜向反射。光粒子在FD坐标系下以天顶角θ和方位角ϕ入射到叶片上，被叶片镜向反射，则C点坐标为：

式中，f表示由坐标系FD到坐标系fL的角度变换，F表示由坐标系fL到坐标系FD的角度变换。

（3）朗伯散射。模型中叶片前/后向散射以及土壤层散射均视为朗伯散射，则C点坐标为

式中，ϕ= 2πδ2，θ'为坐标系fL下光粒子散射天顶角，当发生前向散射时θ'= 0.5πδ1，当发生后向散射时，θ'= π −0.5πδ1，δ1和δ2均为0—1 的均匀随机数。

2.2.3 直射光与散射光

辐射光源为直射光情况时，则设定光粒子从冠层上方不同位置处以固定的天顶角和方位角入射；如果辐射光源为天空散射光，则设定光粒子从冠层上半球空间的不同位置以随机的天顶角和方位角入射到冠层。通过设定两种情形下发射光粒子数量来调整直射光强度和散射光强度的比例。

2.2.4 能量统计

光粒子在各波段上的初始能量均设为单位1，当光粒子进入冠层并与叶片发生碰撞时，设第n次碰撞前光粒子能量为p(n−1)，碰撞中被叶片吸收的能量为pa(n)，碰撞后剩余能量被分为两个部分，其一为光粒子继续传输的能量p(n)，其二为由观测方向冠层空隙逃逸出冠层上边界被观测器直接捕获的能量pu(n)。

冠层顶部的虚拟半球同样被划分为m个等面积的面元，观测器在每个面元处收集光粒子能量，收集到的总能量P为

式中，N为发射总光粒子数量，第一项为光粒子与叶片发生的所有碰撞中Pu的累加值，第二项为沿观测方向逃逸出冠层顶后被观测器捕获的光子总能量。

在(θ，ϕ)观测方向的场景BRF可表示为

3 模型验证与分析

3.1 实验数据

实验观测数据取自2010 年6 月—8 月在中国科学院怀来遥感综合试验站测量的玉米多角度光谱数据，玉米植株平均行间距为58 cm，株间距为38 cm，实验测量了玉米3 个具有代表性的生长期（离散期、垄行期、均一期）的垂直垄方位的多角度BRF 观测数据。观测时段均选择太阳位于垂直垄方位进行，观测天顶角从−60°—60°，间隔10°。整个方位观测时间大约持续20 min，观测期间太阳方位约有5°变化。观测仪器为ASD 光谱仪（400—2500 nm），视场角25°。测量的配套数据包括植株高度、叶展宽和单株叶片总面积等数据，其中，总面积可换算叶面积指数LAI（Leaf Area Index）。实验中利用ASD 光谱仪测量了土壤组分的反射光谱数据，配合使用积分球测量了叶片的反射率和透过率，主要参数见表1和表2。

表1 主要观测参数Table 1 Main observation parameters

表2 组分反射率/透过率数据Table 2 Component reflectance/transmittance data

3.2 模型参数设置

考虑玉米概率株构建的便捷性，如图3（a）所示，在本文中玉米概率株株冠被简化为圆锥与圆台的组合体，圆台下底面半径为其上底面半径的一半，圆锥高度占株高的20%，圆台高度占株高的70%；忽略玉米叶片形状对叶片体密度分布的影响，并设定株冠各水平层内叶片体密度最大值相同而叶片体密度沿径向按余弦规律递减至0。由于叶倾角沿径向和垂直向的分布均存在变化并且变化规律随生长期的不同而不同，在概率株构建时很难详细考虑叶片角度的真实分布。基于对垄行期玉米叶片角度测量数据的分析，本文将概率株的叶倾角分布简化为如图3（b）所示的梯形分布，叶倾角30°—60°出现概率最大且相等，小于30°或大于60°的叶倾角出现概率线性递减至最大概率的25%。考虑光谱仪的视场角，设定收集光粒子能量的面元面积占虚拟半球面积的1/50。模型中玉米秆组分的反射率用叶片反射率代替，透过率为0。按照表2 设定场景组分在红光波段（680 nm）和近红外波段（850 nm）的反射率与透过率参数，模拟光粒子数N设定为2×105，太阳直射光与散射光比例设定为20∶1。图4 展示了观测样地的3个生长期的LAI空间分布。

图3 玉米概率株及叶片体密度分布Fig.3 Schematic diagram of corn probability plant and leaf density distribution

图4 玉米不同生长期LAI空间分布Fig.4 LAI distribution of corn in different growth periods

3.3 模型验证

利用本文模型对离散期、垄行期和均一期的玉米冠层太阳主平面（垂直垄向）的BRF 进行了以5°为间隔的5次模拟，结果如图5所示。由于5次BRF 模拟中概率株坐标在垄向和垂直垄向均有0—2 cm 的随机变化，以及模拟中光粒子总数设定偏小导致BRF的模拟结果有一定的发散，但BRF最大变化量小于0.03，这说明模型具有较好的稳定性。

图5 玉米冠层BRF模拟结果与实测值对比Fig.5 Comparison of corn canopy BRF between simulation results and measured values

冠层叶片对红光有强烈吸收作用，在该波段叶片组分反射率低于土壤组分反射率，随着冠层LAI的增加，视场内叶片组分占比越来越大，导致红光波段地表反射率下降，0°观测角度下地表发射率从0.15（离散期）下降到0.05（均一期）。近红外波段叶片组分反射率高于土壤组分反射率，随着冠层LAI的增加，地表近红外波段反射率整体上升，0°观测角度下地表发射率从0.29（离散期）增加到0.41（均一期）。

3个生长期的冠层BRF 模拟结果均有明显的热点效应，而‘热点’角度下BRF 的实测数值均低于模拟值，红光波段平均低0.023，近红外波段平均低0.026，这主要是由于ASD 光谱仪较大的观测视场角对BRF 随角度的变化具有一定的平滑作用，另外多角度观测装置在热点方向的阴影同样会降低‘热点’BRF的观测值。

在离散期，地表前向BRF 模拟结果相比实测值偏高，60°观测角度下，地表红光波段和近红外波段的BRF 模拟值均比观测值高出0.04 以上，这主要是由于离散期的玉米植株较小，观测视场内以土壤组分为主，而模拟中把土壤介质假设为朗伯体，这与实际不符。随着LAI的增大，观测视场内土壤组分占比下降，在垄行期和均一期BRF 模拟值偏高的现象消失，同时叶片对近红外波段的强散射使冠层BRF的碗边效应更加明显。

玉米冠层3 个生长期5 次BRF 模拟结果随观测角度的拟合曲线与实测数据保持了较好的一致性变化趋势，红光波段模拟值与实测值的RMSE为0.0085，R2为0.96，近红外波段模拟值与实测值的RMSE 为0.013，R2为0.96。验证结果表明，本文模型具备准确模拟不同生长期玉米冠层BRF 的能力。影响模型模拟精度的因素除上述分析中提到的模型光粒子N的设置、观测视场角过大、观测装置阴影以及土壤组分的朗伯散射近似等原因外，概率株模型构建的逼真度以及实验观测视场尺度大小等因素的影响也是不可忽略的。一方面概率株是基于真实植株相关参数的统计规律构建的，是对大量真实植株的一种近似，因此概率株构建中的近似与假设条件的引入以及植株相关参数的实测数据本身的测量误差、数据量不富足或代表性不突出等均会影响概率株构建的逼真度，从而影响冠层BRF 的模拟精度，比如，本文玉米概率株构建中对叶片体密度分布以及叶片角度分布均进行了一定的近似和假设，因此引入的误差是不可忽略的；另一方面，地基BRF 观测实验中观测视场尺度有限，视场内植株数量少，而概率株冠层场景是大场景植株冠层的等效，两者间存在的不确定性差异也会影响模型的验证结果。

3.4 光粒子数量对BRF的影响

在相同的参数设置下，将光粒子数分别设定为2×105和5×106进行均一期冠层BRF 随观测方位角和天顶角变化的模拟，结果如图6所示。当光粒子数为2×105时，BRF模拟用时3.104 min，由于模拟结果有一定的随机性导致BRF分布图有‘斑点’出现；当光粒子数设为5×106时，模拟用时27.531 min，模拟结果的随机性减弱，BRF 分布图的‘斑点’基本消失。由不同光粒子数对冠层BRF 模拟结果可见，随着光粒子数增加，本文模型BRF 模拟结果稳定性有显著提高。

图6 不同光粒子数设置下的玉米冠层BRF分布Fig.6 BRF distribution of corn canopy with different light particles number

3.5 与SAIL模型的对比

冠层二向反射率模型SAIL （Scattering by Arbitrarily Inclined Leaves）被广泛应用于水平均一植被冠层的BRF 模拟中（Jacquemoud，1993），为实现与SAIL 模型的对比，将概率株设定为叶片体密度在概率株平均占地立体空间内均匀分布；利用PROSPECT 模拟的叶片反射率和透过率作为SAIL 模型和本文模型的叶片光谱特性参数。PROSPECT 模型叶片结构参数N设定为1.5，叶绿素含量为40%，胡萝卜素含量为8%；SAIL 模型热点系数0.01，土壤参数为0.5；在太阳天顶角20°，叶片角度分布均为球形分布的条件下，模拟了冠层LAI分别为2.5和0.5时的冠层太阳主平面和垂直主平面BRF 随观测角度的变化，模型对比结果如图7所示。

图7 SAIL模型与本文模型的BRF模拟结果对比Fig.7 Comparison results of simulated BRF between SAIL model and the model developed in this paper

水平均一冠层BRF 的模拟对比结果显示，对不同的LAI 冠层或不同的观测方位，本文模型与SAIL 模型的BRF 模拟值随观测角度的变化均保持了较高的一致性（R2均大于0.83），这说明本文模型在光能粒子传输与能量计算方面的是合理准确的。两者对太阳主平面方位BRF 的模拟结果一致性略低于垂直主平面方位，两者差别主要体现在本文模型模拟的太阳主平面冠层BRF 值略高于SAIL 模型结果，尤其是热点角度至后向最大观测角度范围内本文模型BRF 下降幅度和速度明显小于SAIL 模型，这主要是因为本文模型计算的是一定观测立体角内的平均BRF，对BRF 随角度变化具有一定的平滑作用，而SAIL 模型的BRF 并没有观测立体角的考虑。另外，SAIL 模型对冠层LAI为2.5 的垂直太阳主平面近红外波段BRF 的模拟结果为非对称分布，本文模型的垂直太阳主平面的BRF 明显具有轴对称特征，与均匀冠层BRF 理论分布更加吻合。

4 结 论

本文通过概率株来构建三维冠层，基于蒙特卡洛光线追踪构建了冠层辐射传输模型。本文模型有两个显著的特点：（1）用概率株的方法来构建冠层场景不仅能刻画冠层的非均一分布而且可考虑单植株立体空间内叶片的非均一分布；（2）模型融合了蒙特卡洛技术和冠层孔隙率的计算，基于蒙特卡洛方法模拟双重属性光粒子在冠层内的传输行为，通过冠层孔隙率计算实现对沿观测向冠层空隙逃逸的光能量的单独统计，这对模型计算稳定性的提高起到了关键作用，光粒子数在106量级就能实现稳定的模拟结果。

以不同LAI的水平均一冠层为例，进行了本文模型与SAIL 模型在冠层太阳主平面和垂直主平面BRF 的模拟对比分析，对比结果表明本文模型在光能粒子传输与能量计算方面是准确的，BRF 角度变化模拟结果更加合理。以具有典型垄行种植结构的玉米冠层为例进行了不同生长期冠层多角度BRF 的模拟，并与观测数据进行了比对。在玉米冠层离散期、垄行期和均一期的BRF 模拟结果与实测数据均保持了较好的一致性变化趋势，红光波段模拟值与实测值的RMSE 为0.0085，R2为0.96，近红外波段模拟值与实测值的RMSE 为0.013，R2为0.96。验证结果表明，本文模型具备准确模拟不同生长期玉米冠层BRF 的能力，为玉米冠层BRF 的研究提供了一个有效模拟工具。另外，一般农作物以及阔叶树种等均可看成具有中心对称特征的植株，本模型可适用于该类型植被冠层BRF 的模拟，由概率株和场景参数可确定如生物量、LAI等冠层参数的理论值，因此，本文模型可为植被生物量以及LAI的定量化反演与算法验证提供数据支撑。

本文模型对非均一性冠层场景的构建较真实场景模型简单，与SAIL 模型相比能够模拟植被冠层异质性对BRF的影响，与传统MC类辐射传输模型相比模拟的稳定性更高。本文发展的模型还存在一些不足，需要我们下一步继续改进：（1）概率株模型的构建需要真实植株组分分布规律的先验知识作为支撑，目前概率株模型仅适用于具有统计意义下的中心对称特征的植被类型，而对于非中心对称的植被类型以及有各向异性枝杆的植被类型，如何利用概率株进行该类冠层组分空间分布的有效刻画，还有待进一步研究；（2）模型中关于土壤和叶片组分的二向散射特性进行了朗伯的简化处理，这也是模型需要后续改进的一个方面；（3）本文模型仅使用了玉米冠层的BRF 观测数据进行了初步的验证，未来还需要与其他三维冠层辐射模型进行对比，收集更多植被冠层光谱观测数据对本文模型进行更全面的验证。