CN103413014A - 一种基于二向性反射分布函数(brdf)原型反演地表反照率的方法 - Google Patents

Abstract

本发明涉及一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法。通过应用BRDF原型对多角度观测数据进行循环迭代拟合，求出最小均方根误差(RMSE)，提出基于BRDF原型反演地表反照率的新方法，以黑河地区机载红外广角双模式成像仪(WiDAS)实验数据为例，应用BRDF原型作为二向性反射先验知识，约束大观测天顶角下二向性反射形状的剧烈变化，克服极端反射率变化所引起的反照率的不确定性。该方法不同于目前业务化运行的基于核驱动模型反演地表反照率的反演算法，为空间采样能力不足的传感器提供了一种反演地表反照率的解决方案。本发明在空间信息技术领域，尤其在定量遥感方面具有重要的应用价值。

Description

一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法

一、技术领域

本发明涉及一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法，属于空间信息技术领域。

二、背景技术

地球表面各向异性反射是自然界中物体对电磁波反射的基本宏观现象，地表二向性反射分布函数的研究，是光学定量遥感研究领域重要组成部分。70年代以来，卫星遥感主要采取垂直观测方式，获取地表植被冠层的生物物理参数。随着遥感技术的发展，植被与光辐射之间相互作用机理的深入研究，以及所面临的各种新需求，人们越来越清楚地认识到，在推断地表关键参数方面(如地表反照率)，多角度遥感提供了额外有效的信息源，起到了不可替代的作用。近20年来，随着多角度星载传感器的成功发射，对地研究已积累了大量的多角度观测数据，但受空间采样能力的限制，大多数机载或星载传感器通常不能像地表多角度测量一样，获取足够的空间采样数据，有些稀疏多角度观测数据由于空间信息量的不足，常常会导致反演的失败，因此，如何有效利用这些信息量不足的稀疏多角度观测数据反演地表反照率，是空间信息技术领域亟待需要解决的关键问题之一，也是定量遥感领域各国科学家正在积极研究探索的科学问题之一。

目前，应用稀疏多角度观测数据反演地表反照率，通常依赖于一个二向性反射数据库，早期的二向性反射数据库试图建立地表生态类型和地表二向性反射变化之间的规律，但是，有研究表明地表生态类型二向性反射变化类内方差通常会大于类间方差，这为该方法实际应用带来很大的理论障碍；目前，MODIS采用像元对像元的全球BRDF数据库，该数据库基于高质量历史数据累积，这个数据库本身存在的问题在于：(1)需要假设不同年份相同时期的BRDF变化稳定一致，通过多年数据累积和填补，获得相同时期高质量BRDF数据，而此假设并不能完全成立；(2)即使在该假设成立的情况下，获得全球完整的多角度数据有很大困难，如，在赤道区域和高纬地区，由于云和太阳天顶角的影响，很难获得高质量的BRDF数据；(3)数据量大，即，该方法产生的BRDF数据库，通常有几百个GB且不能获得BRDF变化的基本规律，虽然MODIS传感器依靠这样一个量身定做的数据库获得全球反照率的量反演质量，但这个数据库在用于其它多角度传感器时，由于空间尺度和时间尺度的差异存在很大的困难，如在973项目支持下，对早期黑河区域飞行中获取的机载WiDAS稀疏多角度数据的反演，这个数据库应用很困难。

本发明基于我们早期构建的一个简洁的BRDF函数原型库，开发了一种应用稀疏多角度数据反演地表反照率的新方法，该方法不需要考虑地表生态类型，以及空间尺度和时相尺度的变化，通过应用构建的BRDF原型，对不同波段的稀疏多角度数据进行迭代拟合，获得不同原型和稀疏多角度数据的最小二乘拟合误差(RMSE)，通过选取最小RMSE所对应的原型作为最优BRDF原型，然后，求出该原型所对应的黑天空和白天空反照率，通过野外同步测量的天空散射光比例，计算真实地表反照率，最后，将反演结果与地表反照率表实测反照率进行对比，验证算法精度，结果表明，本发明对于稀疏多角度数据的反演，与不使用BRDF原型约束的算法，其精度有明显的改善，尤其对全反演算法失败的情况改善效果显著，体现了本发明的重要性和实用价值。

三、发明内容

1、目的：本发明目的是利用提取的BRDF原型反演机载WiDAS地表反照率，由于机载WiDAS传感器稀疏采样观测角度窄，导致在利用核驱动模型进行全模型无约束反演时，在大观测角度下缺少约束，使外延的BRDF形状变化剧烈，从而导致反演的地表反照率错误，本发明利用BRDF原型进行约束，通过模型拟合BRDF原型与多角度观测，选择最优的BRDF原型，很好地结合了地表多角度观测与地表BRDF先验知识，成功反演机载WiDAS地表反照率。本发明在定量遥感方面，尤其对利用BRDF先验知识反演地表反照率等关键地表参数方面提供了约束算法，有重要的应用价值。

2、技术方案：本发明涉及一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法，具体步骤如下(如图1)：

步骤一：针对各原型求模型和观测的最小拟合误差

对于稀疏多角度采样数据，BRDF空间采样不足会导致多角度观测信息量不足，因此，用全模型反演地表反照率，由于在大角度观测方向缺少约束，会高估或低估反照率，BRDF原型提供了一种有效约束，通过对稀疏采样数据用BRDF原型进行约束拟合，可综合考虑传感器实际观测和地表各向异性反射变化，拟合过程可表示如下：

假定一组有n个观测的多角度观测数据B：B＝ρ₀(Ω_i，0；Ω_r，0)，ρ₁(Ω_i，1；Ω_r，1)...ρ_n-1(Ω_i，n-1；Ω_r，-1)，若与其对应的BRDF原型为BRDF′，对应相同观测位置的原型反射率数据表示为B′：B′＝ρ′₀(Ω_i，0；Ω_r，0)，ρ′₁（Ω_i，1；Ω_r，1)...ρ′_n-1(Ω_i，n-1；Ω_r，n-1)，取调整系数a，通过平移BRDF原型，使得B和aB′的差异最小，采用最小二乘拟合法最小化拟合误差e²

$e^2=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{({\mathrm{\ρ}}_i-a{{\mathrm{\ρ}}^{\′}}_i)}^2---\left(1\right)$

求调整系数a

$a=\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_j\×{{\mathrm{\ρ}}^{\′}}_j/\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\left({{\mathrm{\ρ}}^{\′}}_j\right)}^2---\left(2\right)$

这样，数据B的BRDF可表示为

BRDF＝a×BRDF′         (3)

拟合误差RMSE可以表示为

$\mathrm{RMSE}=\sqrt{\frac{1}{n-1}\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{({\mathrm{\ρ}}_k-a{{\mathrm{\ρ}}^{\′}}_k)}^2}---\left(4\right)$

步骤二：选取最优原型求其光谱固有反照率

用6个原型分别拟合观测数据，计算各原型的拟合误差RMSE，选择拟合观测数据RMSE最小的原型作为数据的先验BRDF，求该原型的方向-半球和半球-半球积分如下：

$h_{\mathrm{BRDF}\′}\left({\mathrm{\θ}}_s\right)=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}{\∫}_0^{\mathrm{\π}/2}\mathrm{BRDF}\′({\mathrm{\θ}}_v,{\mathrm{\θ}}_s,\mathrm{\Δ\φ})\sin \left({\mathrm{\θ}}_v\right)\cos \left({\mathrm{\θ}}_v\right)d{\mathrm{\θ}}_v\mathrm{d\φ}---\left(5\right)$

对h_BRDF′(θ_s)在太阳入照半球进一步积分，得半球-半球积分

$H_{\mathrm{BRDF}\′}(\mathrm{\θ},\mathrm{\λ})=2{\∫}_0^{\mathrm{\π}/2}{h}_{\mathrm{BRDF}\′}\left({\mathrm{\θ}}_s\right)\sin \left({\mathrm{\θ}}_s\right)\cos \left({\mathrm{\θ}}_s\right)d{\mathrm{\θ}}_s---\left(6\right)$

通过式(3)确定调整系数a，则黑天空反照率(Black Sky Albedo，BSA)和白天空反照率(White SkyAlbedo，WSA)与选定原型的上述积分值之间的关系可以表示为：

BSA(θ，λ)＝α×h_BRDF′(θ，λ)      (7)

WSA(λ)＝a×H_BRDF′(λ)      (8)

步骤三：产生地表宽波段真实反照率

根据已获得的光谱固有反照率，即白天空反照率(WSA)和黑天空反照率(BSA)，考虑到大气状态，通过测量天空散射光和直接光的比值，获得光谱真实反照率

α_blue(θ_i，λ)＝SKYL(θ_i，λ)×WSA(λ)+(1-SKYL(θ_i，λ))×BSA(θ_i，λ)    (9)

通过建立WiDAS传感器波段与一个标准波段(如MODIS前7个波段)之间的关系，将WiDAS传感器的光谱反照率转换成标准波段的光谱返归率，通过窄波段到宽波段的转换公式，获得地表宽波段的真实反照率

α_short＝0.160α₁+0.291α₂+0.243α₃+0.116α₄

(10)

+0.112α₅+0.081α₇-0.0015

最后，通过地面验证和全模型无约束反演进行比较，验证算法精度。

3、优点及功效：本发明开发了一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法，通过BRDF原型约束核驱动模型，对稀疏采样多角度观测数据实施了基于先验知识约束的反照率反演。本发明方法适应于传感器稀疏采样的观测数据集，逻辑清楚，适应性强，综合考虑了地表实际观测值和地表BRDF先验知识积累，在区域和全球尺度上，针对多角度采样能力不足，或由于云的影响所造成多角度观测信息量不足的情况，快速准确地反演地表反照率等地表关键参量有理论和应用价值。

四、附图说明

图1本发明方法的流程示意图

图2(a)黑河盈科绿洲实验区机载WiDAS无约束全反演反照率(2008年6月29号)

图2(b)黑河盈科绿洲实验区机载WiDAS原型约束反演反照率(2008年6月29号)

图3(a)无约束全反演的极小值反照率所对应的角度采样和主平面BRDF形状(近红外)

图3(b)无约束全反演的合理值反照率所对应的角度采样和主平面BRDF形状(近红外)

图3(c)无约束全反演的极大值反照率所对应的角度采样和主平面BRDF形状(近红外)

图4验证点地表实时通量塔和反照率表分布图

图5两种方法(无约束全反演和约束原型反演)产生反照率的验证

五、具体实施方式

本发明涉及一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法，该方法具体步骤如下(如图1)：

步骤一：针对各原型求模型和观测的最小拟合误差

对于稀疏多角度采样数据，BRDF空间采样不足会导致多角度观测信息量不足，因此，用全模型反演地表反照率，由于在大角度观测方向缺少约束，会高估或低估反照率，BRDF原型提供了一种有效约束，通过对稀疏采样数据用BRDF原型进行约束拟合，可综合考虑传感器实际观测和地表各向异性反射变化，拟合过程可表示如下：

假定一组有n个观测的多角度观测数据B：B＝ρ₀(Ω_i，0;Ω_r，0)，ρ₁(Ω_i，1；Ω_r，1)...ρ_n-1(Ω_i，n-1；Ω_r，n-1)，若与其对应的BRDF原型为BRDF′，对应相同观测位置的原型反射率数据表示为B′：B′=ρ′₀(Ω_i，0；Ω_r，0)，ρ′₁(Ω_i，1；Ω_r，1)...ρ′_n-1(Ω_i，n-1；Ω_r，n-1)，取调整系数a，通过平移BRDF原型，使得B和aB′的差异最小，采用最小二乘拟合法最小化拟合误差e²

$e^2=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{({\mathrm{\ρ}}_i-a{{\mathrm{\ρ}}^{\′}}_i)}^2---\left(1\right)$

求调整系数a

$a=\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_j\×{{\mathrm{\ρ}}^{\′}}_j/\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\left({{\mathrm{\ρ}}^{\′}}_j\right)}^2---\left(2\right)$

这样，数据B的BRDF可表示为

BRDF＝a×BRDF′       (3)

拟合误差RMSE可以表示为

$\mathrm{RMSE}=\sqrt{\frac{1}{n-1}\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{({\mathrm{\ρ}}_k-a{{\mathrm{\ρ}}^{\′}}_k)}^2}---\left(4\right)$

步骤二：选取最优原型求其光谱固有反照率

用6个原型分别拟合观测数据，计算各原型的拟合误差RMSE，选择拟合观测数据RMSE最小的原型作为数据的先验BRDF，求该原型的方向-半球和半球-半球积分如下：

$h_{\mathrm{BRDF}\′}\left({\mathrm{\θ}}_s\right)=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}{\∫}_0^{\mathrm{\π}/2}\mathrm{BRDF}\′({\mathrm{\θ}}_v,{\mathrm{\θ}}_s,\mathrm{\Δ\φ})\sin \left({\mathrm{\θ}}_v\right)\cos \left({\mathrm{\θ}}_v\right)d{\mathrm{\θ}}_v\mathrm{d\φ}---\left(5\right)$

对h_BRDF′(θ_s)在太阳入照半球进一步积分，得半球-半球积分

$H_{\mathrm{BRDF}\′}(\mathrm{\θ},\mathrm{\λ})=2{\∫}_0^{\mathrm{\π}/2}{h}_{\mathrm{BRDF}\′}\left({\mathrm{\θ}}_s\right)\sin \left({\mathrm{\θ}}_s\right)\cos \left({\mathrm{\θ}}_s\right)d{\mathrm{\θ}}_s---\left(6\right)$

通过式(3)确定调整系数a，则黑天空反照率(Black Sky Albedo，BSA)和白天空反照率(White Sky Albedo，WSA)与选定原型的上述积分值之间的关系可以表示为：

BSA(θ，λ)＝a×h_BRDF′(θ，λ)      (7)

WSA(λ)＝a×H_BRDF′(λ)      (8)

步骤三：产生地表宽波段真实反照率

根据已获得的光谱固有反照率，即白天空反照率(WSA)和黑天空反照率(BSA)，考虑到大气状态，通过测量天空散射光和直接光的比值，获得光谱真实反照率

α_blue(θ_i，λ)＝SKYL(θ_i，λ)×WSA(λ)+(1-SKYL(θ_i，λ))×BSA(θ_i，λ)      (9)

通过建立WiDAS传感器波段与一个标准波段(如MODIS前7个波段)之间的关系，将WiDAS传感器的光谱反照率转换成标准波段的光谱返归率，通过窄波段到宽波段的转换公式，获得地表宽波段的真实反照率

α_short＝0.160α₁+0.291α₂+0.243α₃+0.116α₄

(10)

+0.112α₅+0.081α₇-0.0015

最后，通过地面验证和全模型无约束反演进行比较，验证算法精度。

实施例1：

在一台配置有Intel(R)Core^TM2.5GHz4处理器，4G内存，ATI Radeon HD5670图形显卡的计算机上进行了实施，以机载WiDAS多角度观测数据为例，采用本方法，进行稀疏采样多角度数据的反照率反演(图1)。

对比无约束全模型反演结果和有约束原型反演结果(图2(a)和图2(b))，可以看出无约束全模型反演结果(图2(a))明显存在大范围小于0的反照率值，这是由于WiDAS传感器采样范围小，在大观测角度下由于缺少约束，使大角度外延的反射率明显小于0(图3(a)中点划线)造成的，BRDF原型的约束很好克服了这一点(图3(a)中实线)，对于无约束全模型反演结果(图2(a))存在少量大于1的反照率值，在图3(c)中也给出类似的解释，当然，无约束全模型反演结果也存在反演结果合理的值，这些值所对应的BRDF形状和所用的BRDF原型很接近，如(图3(b))。

进一步用地表实时布设的反照率表(图4)进行精度验证(图5)，可以看出，无约束全模型反演结果不确定性很大，反照率值的变化范围约在-0.1～0.15，对同一种地表类型，这种变化显然是错误的，但有约束原型反演结果很稳定，变化范围在0.1～0.15，反演的反照率值绝对精度在0.02以内，很好地满足反照率精度要求(图5)。

Claims (1)

1.一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法，步骤包括：

步骤一：针对各原型求模型和观测的最小拟合误差

对于稀疏多角度采样数据，BRDF空间采样不足会导致多角度观测信息量不足，因此，用全模型反演地表反照率，由于在大角度观测方向缺少约束，会高估或低估反照率，BRDF原型提供了一种有效约束，通过对稀疏采样数据用BRDF原型进行约束拟合，可综合考虑传感器实际观测和地表各向异性反射变化，拟合过程可表示如下：

假定一组有n个观测的多角度观测数据B：B＝ρ₀(Ω_i，0；Ω_r，0)，ρ₁(Ω_i，1;Ω_r，1)...ρ_n-1(Ω_i，n-1；Ω_r，n-1)，若与其对应的BRDF原型为BRDF′，对应相同观测位置的原型反射率数据表示为B′：B′=ρ′₀(Ω_i，0；Ω_r，0)，ρ′₁(Ω_i，；Ω_r，1)...ρ′_n-1(Ω_i，n-1；Ω_r，n-1)，取调整系数a，通过平移BRDF原型，使得B和aB′的差异最小，采用最小二乘拟合法最小化拟合误差e²

$e^2=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{({\mathrm{\ρ}}_i-a{{\mathrm{\ρ}}^{\′}}_i)}^2---\left(1\right)$

求调整系数a

$a=\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_j\×{{\mathrm{\ρ}}^{\′}}_j/\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\left({{\mathrm{\ρ}}^{\′}}_j\right)}^2---\left(2\right)$

这样，数据B的BRDF可表示为

BRDF＝a×BRDF′       (3)

拟合误差RMSE可以表示为

$\mathrm{RMSE}=\sqrt{\frac{1}{n-1}\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{({\mathrm{\ρ}}_k-a{{\mathrm{\ρ}}^{\′}}_k)}^2}---\left(4\right)$

步骤二：选取最优原型求其光谱固有反照率

用6个原型分别拟合观测数据，计算各原型的拟合误差RMSE，选择拟合观测数据RMSE最小的原型作为数据的先验BRDF，求该原型的方向-半球和半球-半球积分如下：

$h_{\mathrm{BRDF}\′}\left({\mathrm{\θ}}_s\right)=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}{\∫}_0^{\mathrm{\π}/2}\mathrm{BRDF}\′({\mathrm{\θ}}_v,{\mathrm{\θ}}_s,\mathrm{\Δ\φ})\sin \left({\mathrm{\θ}}_v\right)\cos \left({\mathrm{\θ}}_v\right)d{\mathrm{\θ}}_v\mathrm{d\φ}---\left(5\right)$

对h_BRDF′(θ_s)在太阳入照半球进一步积分，得半球-半球积分

$H_{\mathrm{BRDF}\′}(\mathrm{\θ},\mathrm{\λ})=2{\∫}_0^{\mathrm{\π}/2}{h}_{\mathrm{BRDF}\′}\left({\mathrm{\θ}}_s\right)\sin \left({\mathrm{\θ}}_s\right)\cos \left({\mathrm{\θ}}_s\right)d{\mathrm{\θ}}_s---\left(6\right)$

通过式(3)确定调整系数a，则黑天空反照率(Black Sky Albedo，BSA)和白天空反照率(White SkyAlbedo，WSA)与选定原型的上述积分值之间的关系可以表示为：

BSA(θ，λ)＝a×h_BRDF′(θ，λ)      (7)

WSA(λ)＝a×H_BRDF′(λ)      (8)

步骤三：产生地表宽波段真实反照率

根据已获得的光谱固有反照率，即白天空反照率(WSA)和黑天空反照率(BSA)，考虑到大气状态，通过测量天空散射光和直接光的比值，获得光谱真实反照率

α_blue(θ_i，λ)＝SKYL(θ_i，λ)×WSA(λ)+(1-SKYL(θ_i，λ))×BSA(θ_i，λ)    (9)

通过建立WiDAS传感器波段与一个标准波段(如MODIS前7个波段)之间的关系，将WiDAS传感器的光谱反照率转换成标准波段的光谱返归率，通过窄波段到宽波段的转换公式，获得地表宽波段的真实反照率

α_short＝0.160α₁+0.291α₂+0.243α₃+0.116α₄

(10)

+0.112α₅+0.081α₇-0.0015

最后，通过地面验证和全模型无约束反演进行比较，验证算法精度。

Images