一种计算密实性颗粒介质光散射特性的方法

摘要

本发明提供一种计算密实性颗粒介质光散射特性的方法，将密实性介质的光散射特性——二向反射分布函数(BRDF)分为两部分来计算，即面散射BRDF和体散射BRDF。通过波解析理论计算得到散射介质面散射、通过辐射传输理论计算得到体散射的二向反射分布特性BRDF，两部分计算之和为密实性颗粒介质光散射特性。本发明首次将辐射传输理论与波解析理论相结合，相对于以前的密实性介质光散射模型而言，充分考虑了介质表面形状对光散射特性的作用，增加了刻画散射介质表面粗糙程度的参数，提高了计算精度。

权利要求

1.一种计 算密实性 颗粒介 质光散射 特性的方 法，介质 光散射特 性其中，θi为入射光天顶角，为入射光方位角，θs为观测天顶角，为观测方位角，λ为入射光波长； 为介质面散射的二向反射分布特性， 为介质体散射的二向反射分布特性；

所述介质面散射的二向反射分布特性其中，

k＝2π/λ，π为圆周率，T为介质表面的特征长度，H＝3-D，D为介质表面的分维数，0＜H＜1Γ()为伽马函数，

vz＝k(cosθs+cosθi)，或

或

或

Rh0为水平极化菲涅耳反射系数，Rv0为垂直极化菲涅耳反射系数；

所述介质体散射的二向反射分布特性其中，

ω为散射介质单次散射反照率，μ0＝cosθi，u＝cosθs，bn为相函数勒让德展开系数，Pn()为勒让德多项式，g为μ0或u；

2.根据权利要求1所述的计算密实性颗粒介质光散射特性的方法，其特征在于，在θs＝0°时考虑遮蔽效应，则 介质面散射的二向反射分布特性tan(θ0)表示粗糙面的均方根斜率，erfc()为余误差函数。

说明书

一种计算密实性颗粒介质光散射特性的方法

技术领域

[0001] 本发明涉及一种计算密实性介质光散射特性的方法，计算的光谱范围为可见光波段即0.4微米到2.5微米。

背景技术

[0002] 介质的光散射特性是指光波照射到介质上后介质对光波的影响(包括光的反射和折射以及光被吸收)。对于密实性介质，可以通过研究介质的二向反射特性(具体的来说是介质的二向反射分布函数BRDF)来研究介质光散射特性。影响介质二向反射特性的主要因素有介质的组成成分、各个组分的内部晶体结构、介质的表面形状。因此，研究介质的二向反射特性可以达到两个目的：一方面，在已知物质结构、属性、形状的特性的情况下，得到介质的二向反射光谱曲线；另一方面，由遥感或实验室测量得到的部分光谱数据通过反演得到二向反射模型的参数，根据这些参数我们可以对介质的组分等进行定量的分析。另外，将参数带入模型中，根据模型可得到介质的整个光谱曲线。

[0003] 随着高光谱、航天及遥感定量反演技术的发展，物体的二向性反射特性不应被忽视。地物二向性反射模型的建立与完善不仅可以仿真地物实测光谱，而且也是定量反演地表参数的重要理论基础。根据地物的光谱特征可以直接识别地物组成以至地物成分，光谱重建是高光谱识别和量化反演的前提和基础。因此，发展二向性反射模型对于提高地物成分反演、以及光谱解混的精度意义重大。

[0004] 为计算密实性散射介质(行星表面风化层、裸露地表、岩矿粉末等)的光学散射特性，在钱德拉萨卡(Chandrasekhar)辐射传输方程基础上，研究者建立了不同类型的基于辐射传输理论的二向性反射模型，目前最常用的是哈帕克(Hapke)模型。

[0005] Hapke模型的发展主要有三个阶段，1981年提出的Hapke模型在计算密实介质对光散射时，假设介质是各向同性的。而自然界中，介质不可能是各向同性的。1993年的Hapke模型引入了粗糙度阴影函数，考虑了表面粗糙度的遮蔽效应；2002的Hapke模型在计算介质散射时，考虑到了各向异性对光散射的影响。得到了更精确的各向异性的光散射的表达式，同时增加了干涉后向效应，并考虑了介质孔隙度对模型的影响。Hapke模型虽然只是辐射传输方程的二流近似解析解，但由于其求解简单、精度基本达到应用需求，并且有着易于反演的优势，因此被广泛应用于行星遥感领域和裸露地表的遥感模型。

[0006] Hapke模型是目前描述密实性散射介质二向性反射的主流模型，但在实际应用中，发现该模型存在一定的缺陷。具体而言，采用Hapke模型模拟密实性散射介质二向性反射光谱以及用模型反演散射介质的光学特性时，发现单次散射反照率大的散射介质模型精度较低，而单次散射反照率低的散射介质模型精度较高；另外，正如前面所说的，影响介质二向反射特性的主要因素有介质的组成成分、各个组分的内部晶体结构、介质的表面形状。而在Hapke模型中没有考虑到密实性介质表面形状对散射光谱的影响即没有考虑到面散射。

发明内容

[0007] 本发明的目的是提供一种计算密实性颗粒介质光散射特性的方法，计算精度高，同时精确反映介质表面的粗糙度。

[0008] 一种计算密实性颗粒介质光散射特性的方法，介质光散射特性其中，θi为入射光天顶角，为入射光方位角，θs为观测天顶角，为观测方位角，λ为入射光波长； 为介质面散射的二向反射分布特性， 为介质体
散射的二向反射分布特性。

[0009] 与Hapke模型相比，本发明的技术效果体现在两个方面：一方面，本发明面散射部分充分的考虑了表面形状对光散射特性的作用，增加了刻画散射介质表面粗糙程度的参数，能够反映实际散射体表面粗糙度的大小。而现有的Hapke模型没有考虑到介质表面形状对光散射特性的影响；另一方面，本发明将辐射传输理论与波解析理论相结合，兼顾了两者的优点，实验仿真表明本发明比Hapke更加准确。

附图说明

[0010] 图1是光与密实性颗粒介质相互作用的剖面示意图。

[0011] 图2是本发明步骤流程图。

[0012] 图3是本发明实施例中散射介质表面几何示意图。图中入射角与x轴夹角为0，所以

[0013] 图4是我们的模型、Hapke模型与离散坐标法(DISORT)数值计算结果的比较图。

[0014] 图5是改变模型输入参数，图4所述各种模型计算结果比较的另一示例图。

具体实施方式

[0015] 本发明提供一种用于计算密实性散射介质的光学散射特性模型，模型的输出为密实性散射介质的二向反射分布函数(BRDF)。如图1所示，介质的光散射包括面散射和体散射。面散射即介质表面对光的作用，体散射是介质内部颗粒对光的作用。因此，计算二向反射分布函数BRDF的步骤包含两部分：首先计算散射介质面散射的BRDF，再计算散射介质体散射的BRDF，两次计算之和为模型输出的BRDF。在面散射中充分考虑了介质表面形状对散射的影响。

[0016] 下面参照附图2对本发明的算法流程详细说明。说明主要分为5部分，1、对BRDF进行简单扼要的介绍；2、简介fBM分形3、介绍单次散射(面散射)的BRDF计算方法；4、介绍多次散射(体散射)的BRDF计算方法。5、仿真比较，将我们的模型及Hapke模型与Disort数值计算进行比较。

[0017] 1、BRDF简介

[0018] 如图2所示，本实施例计算二向反射分布函数(BRDF)的步骤包括：首先计算散射介质面散射的BRDF，再计算散射介质体散射的BRDF，两次计算的和为模型的总的BRDF。

[0019] 二向反射分布函数BRDF，表示观察方向的辐射亮度 与入射方向的辐射照度 的比值。定义式如下：

[0020]

[0021] 式中：θi为入射光天顶角，是入射光方位角，θs是观测天顶角，是观测方位角，λ是入射波长。

[0022] 散射亮度 表示沿着观察方向单位面积、单位立体角上的辐射通量；入射照度 表示入射到介质单位面积上的辐射通量。由电磁场散射理论可得到，[0023]

[0024]

[0025] 式中：ψi、ψs分别表示入射通量和散射通量，A为入射光照射面积；R为散射介质i s表面与光学传感器间的距离；Ep、Epq 分别为入射场和散射场；p、q表示入射波和散射波的极化方式，p、q可以取h或v，其中h表示水平极化，v表示垂直极化；<>表示统计平均值，上标*表示共轭，d表示取单位量，Ω表示立体角。

[0026] 将上述关系式带入BRDF定义式中可以得到式(1)，如下

[0027]

[0028] 2、fBM(布朗)分形模型

[0029] 研究密实性介质的光谱模型时，有两个方面对光谱模型起着至关重要的作用，一方面是介质的本身材质属性，另一方面是介质表面形状。因此我们需要对介质表面进行刻画，换句话就是要对表面进行建模。只有在表面模型确定的情况下，我们才能更好的运用各种方法得到介质的光谱模型。

[0030] 一般情况下，我们使用经典的高斯模型、指数模型等来刻画介质表面，但近年来，国内外的工作者们通过大量的实验研究发现，经典模型并不能很好的描述自然介质表面，而fBM分形能更好的反应真实的自然表面。因此我们使用fBM分形来刻画自然界中密实性介质的表面。

[0031] 用二维fBM来模拟密实性散射介质粗糙面z(x，y)粗糙面时，粗糙面在任意两点2 2H
(x′，y′)，(x″，y″)处的高度差服从期望为0方差为sτ 的高斯分布，表达式如下：

[0032]

[0033] 其中：Q{[z(x′，y′)-z(x″，y″)]＜ζ}表示两点的高度差小于某一给定数值ζ的概率， 表示两点的距离。对于fBM粗糙面，H＝3-D，D为介质表面的分维数。实参量s是与粗糙面的特征长度T有关的物理量，特征长度T为粗糙面上(1-H)
斜率均方根为1的任意两点间的距离，s与T之间满足关系：s＝T 。

[0034] fBM粗糙面具有自仿射性质，即如果令Δz(τ)＝z(x′，y′)-z(x″，y″)，则对H任意的a＞0，Δz(aτ)与aΔz(τ)有相同的统计特征。式(2)是从统计的意义上来刻画fBM分形表面。

[0035] 3、面散射的BRDF的计算方法

[0036] 在计算散射介质面散射的BRDF时，本发明采用fBM(即布朗)分形函数模拟粗糙的散射介质表面，再通过基尔霍夫近似法计算得到面散射的BRDF。

[0037] 图2是本发明实施例中密实性颗粒介质表面几何示意图，从方向 来的平行光入射到介质表面上，利用基尔霍夫近似可得到远场区的散射场表达式：

[0038]

[0039] 式中，入射波波数k＝2π/λ； exp()表示指数函数；π为圆周率；θi是入射天顶角；是入射方位角；θs是观察天顶角；是观察方位角；为入射方向单位矢量；i
为散射方向单位矢量；为入射光照射面积S上的点的矢量坐标；d表示微分；|Ep|为入射场强度； 是与入射方向和散射方向有关的物理量，根据入射波和散射波的极化形式的不同分为四种情况，具体表达形式如下：

[0040] a、入射波和散射波全为水平极化

[0041]

[0042] b、入射波水平极化，散射波垂直极化

[0043]

[0044] c、入射波垂直极化，散射波水平极化

[0045]

[0046] d、入射波和散射波都为垂直极化

[0047]

[0048] Rh0为水平极化菲涅耳反射系数，Rv0为垂直极化菲涅耳反射系数。

[0049] 由(3)式可以计算出平均的散射光强

[0050]

[0051]

[0052] 式中：

[0053]

[0054]

[0055] 称作特征函数，对于哈斯特(Hurst)指数0＜H＜1的fBM分形粗糙面可以推导出，

[0056]

[0057] 将(9)式代入(8)式，并利用坐标变换 可得：

[0058]

[0059]

[0060] 其中 J0(t)为零阶贝赛尔函数，|g|表示求模运算，具体形式如下：

[0061]

[0062] 对(10)式进行积分得到，

[0063]

[0064] Γ()为伽马(gamma)函数。在数值仿真中，需要对无穷求和项进行截断，一般取0＜N≤50。

[0065] 将式(12)代入BRDF的计算式(1)可以得到单次散射的BRDF：

[0066]

[0067] 当考虑遮蔽效应时，我们需要在单次散射的BRDF后乘以分形遮蔽函数。国外学者迈克尔克谢泼德(Micheal K.Shepard)通过对自己的实验数据进行拟合后，给出垂直观测即θs＝0°时，分形遮蔽函数随入射天顶角θi变化的表达式：

[0068]

[0069] 其中，tan(θ0)表示粗糙面的均方根斜率，erfc()为余误差函数。

[0070] 但是该表达式并不能与学者汉诺帕威艾伦(Hannu Parviainen)的实验数据相匹配，我们通过MATLAB2008a软件对Micheal K.Shepard和Hannu Parviainen及自己的数据进行拟合后，得到垂直观测(θs＝0°)时，分形遮蔽函数随入射天顶角θi变化的表达式：

[0071]

[0072] 因此，考虑遮蔽效应时，在垂直观测方向上的分形表面的面散射BRDF的表达为：

[0073]

[0074]

[0075] 4、体散射BRDF的计算方法

[0076] 在钱德拉萨卡(Chandrasekhar)辐射传输方程基础上，并参考2002年的Hapke光散射模型得人计算方法，我们给出了介质体散射的BRDF。具体表达式如下：

[0077]

[0078] 式中：ω为散射介质单次散射反照率；μ0＝cosθi，u＝cosθs；M(μ0，u)的表达式如下：

[0079]

[0080] 式中：

[0081]

[0082]

[0083]

[0084] P(μ0)，P(μ)，可由下式计算

[0085]

[0086]

[0087]

[0088] 式中，An计算式为：

[0089]

[0090] bn为相函数勒让德展开系数，Pn()为勒让德多项式。一般情况下，我们采用非对称的双Henyey-Greenstein相函数q(g)，表达形式如下：

[0091]

[0092] 其中，l是入射方向与观察方向的夹角，c、ξ1、ξ2是表达中给定三个参数，与介质的本身物理属性等有关。

[0093] 将(13)式与(15)求和，我们便得到了密实性散射介质的总的BRDF，表达式如下：

[0094]

[0095]

[0096] 同样，我们可以得到在垂直观测方向上考虑遮蔽效应时，总的BRDF表达如下：

[0097]

[0098]

[0099] 5、仿真比较

[0100] 图4及图5是用我们的模型和Hapke模型计算密实性颗粒介质的BRDF，并与DISORT数值计算结果的比较图。

[0101] DISORT数值计算是求解辐射传输方程的离散坐标法软件模块，被广泛应用于陆地、海洋、大气等与辐射传输相关的领域。作为辐射传输方程的精确解，DISORT模型已经成为验证其它模型准确性的参考标准。

[0102] 图4为本发明模型、Hapke模型及Disort数值计算随θs的变化曲线。我们的模型2
输入参数为：λ＝1um，H＝0.85，ω＝0.94，s ＝0.2，θi＝60°， ε1＝
0.6，ε2＝0.6，c＝0。Hapke模型的输入参数为，θi＝60°， ω＝0.94，ε1＝0.6，ε2＝0.6，c＝0。从图3可以看出改进的Hapke模型计算结果更接近于DISORT数值计算。

[0103] 图5同图4一样，本发明的模型输入参数为：λ＝1um，H＝0.75，ω＝0.94，s2＝0.2，θi＝60°， ε1＝0.75，ε2＝0.75，c＝0。Hapke模型的输入参数为，θi＝60°， ω＝0.94，ε1＝0.75，ε2＝0.75，c＝0。图4也表明改进的Hapke模型比Hapke模型的计算精度要高。