陆表二向反射（BRDF）反演方法研究进展

韩源，闻建光，肖青，鲍云飞，陈曦，刘强，贺敏

1.中国科学院空天信息创新研究院 遥感科学国家重点实验室,北京 100101;2.中国科学院大学 资源与环境学院,北京 100049;3.北京空间机电研究所,北京 100083;4.北京师范大学 全球变化与地球系统科学研究院,北京 100875

1 引言

由于地球表面的反射具有明显的非朗伯特性，当电磁波作用在地球表面时，其地物反射特性在各个方向不同，呈现方向性。为了定量刻画此特性，Nicodemus 等（1977）提出了二向反射分布函 数BRDF（Bidirectional Reflectance Distribution Function）的概念。地表BRDF 的研究对于定量遥感的发展有重要意义，因为其从角度维度上增加了卫星数据表征地物特性的信息量，尤其能表征地表地物的三维结构特征，因此可基于该特征反演地表参数，如地表反照率（Schaaf 等，2002，2008；Muller 等，2012）、叶面积指数（Hasegawa等，2010；Pocewicz 等，2007）、聚集指数（Chen等，2005；Wei 和Fang，2016）等关键陆表参数。另一方面，地表的二向反射特性还可以用于归一化不同角度下得到的反射率，消除因入射角度或者观测角度不同导致观测得到的反射率的差异，这对于归一化植被指数估算（León-Tavares 等，2017）、遥感地表覆盖分类（Guan 等，2020）、长时间序列分析等具有重要意义。

根据文献统计学方法，截止2021 年5 月，使用主题：（“BRDF”OR“Bidirectional Reflectance Distribution Function”）和主题：（“model”OR“inversion”OR“retri-eval”）在Web of Science 数据库中进行检索，结果发现自2010年以来，与BRDF模型或者反演方法有关联的研究文献有1013 篇，内容主要和反照率albedo 反演、叶面积指数LAI反演、遥感产品真实性检验等密切相关。其中，研究的领域主要集中于中低分辨率卫星传感器，如MODIS 传感器、POLDER 传感器，此外，发现对BRDF的研究开始聚焦于山地复杂地表（图1）。对这1013 篇文献统计发现，和BRDF 模型或反演方法相关的论文出版数量呈较平稳的变化趋势，而被引频次呈逐年增加趋势（图2）。

图1 2012年—2021年涉及BRDF模型或反演方法的文献的领域来源Fig.1 The statistics of sources of publications and citations correlated with BRDF model or inversion method from 2012 to 2021

图2 2012年—2021年涉及BRDF模型或反演方法的文献出版和引用情况Fig.2 The statistics of publications and citations correlated with BRDF model or inversion method from 2012 to 2021

图3 BRDF反演方法Fig.3 BRDF inversion methods

准确刻画地表的二向反射特征，离不开高精度的BRDF模型及稳定的反演方法，更需要多角度观测数据的支持。目前，具备多角度观测特征的卫星数据已对外发布了不同尺度的BRDF产品，如MODIS 传感器的BRDF 产品（Lucht 等，2000；Schaaf 等，2002）、POLDER 传感器的BRDF 产品（Lacaze，2006；Lacaze 和Maignan，2006）以及MSG-SEVIRI 传感器的BRDF 产品（van Leeuwen 和Roujean，2002）等。从发布的BRDF 产品算法来看，主要采用了半经验的BRDF 模型，如MOIDS及POLDER BRDF 产品采用了核驱动模型（Lucht等，2000；Schaaf 等，2002）、MFG-MVIRI BRDF产品采用了PRV 模型。围绕BRDF 模型反演中的观测数据不足问题，一些学者提出了扩充观测信息源的方法，如综合多源遥感数据或引入先验知识等（Li 等，2001a；Quaife 和Lewis，2010；Wen等，2017）；另外一些学者提出了控制信息量的方法，即让有限的观测信息尽量分配到目标参数反演中（李小文 等，1997；高峰 等，1998）。此外，针对BRDF模型反演中的观测异常值和观测噪声问题，学者们也提出了不同的解决方法（Zeng 等，2016；Gao等，2002）。对于陆表二向反射（BRDF）而言，发展一套准确稳定的反演方法，对于推动定量遥感的广泛应用具有重要的意义。

本文从BRDF 反演原理出发，指出了常见的BRDF 反演问题，之后在对现有BRDF 反演方法总结概述的基础上，重点介绍了地表BRDF反演的原理、特点和应用实例，并讨论了未来提高陆表BRDF 反演精度的方向。文章除了可以为BRDF 的反演提供方法理论，对以线性模型为基础的其他地学参数的反演同样有借鉴意义，能为今后相关研究提供参考。

2 BRDF反演基本原理

2.1 BRDF遥感反演基本原理

地球表面复杂且范围宽广，所以用遥感手段来反演地表BRDF是最有效的方法。通过模型可以将遥感信号和地表BRDF 联系在一起，BRDF 本质是二项分布函数，该函数可以用核驱动等模型来描述，本文统一将描述BRDF 函数的模型表述为“BRDF 模型”。卫星观测数据、模型以及模型输入参数的关系如式（1）：

所有模型都包含两个空间：数据空间D和模型空间M。对于BRDF 模型，数据空间D由卫星观测的二向反射率因子dBRF（Bidirectional Reflection Factor）组成；模型空间M由模型参数i，v，λ，s组成，其中i，v表示入射和观测方向，λ表示波段，s表示空间结构参数集合。εx表示误差项，包含模型自身误差和观测误差。

从模型空间M到数据空间D，即输入模型参数通过模型模拟观测数据的过程，称为正演过程。反之，通过观测数据来反推模型参数就称为反演过程。在BRDF 反演过程中，首先需要反演得到BRDF 模型的输入参数，然后正演出半球其他方向的BRDF。如果BRDF 模型可表达成简洁的解析形式，如核驱动模型，那么可以直接求解其反模型f(i，v，λ，s)-1，从而得到核系数的解析解。但对于形式十分复杂的物理模型和计算机模型，很难直接推导其反模型，这时候需要采取优化迭代或查找表等方法反演模型参数，即通过调整模型参数正向模拟观测数据，将之与真实观测数据比较，将模拟结果和观测结果最为接近时的模型参数视为最优反演结果。

2.2 BRDF模型介绍

对地表二向反射特性的观测可以追随到20 世纪60 年代（Salomonson，1968），人们发现地表二向反射函数BRDF 能很好的反应地表特征。目前，已发展的BRDF模型可分为物理模型、经验及半经验模型、计算机模拟模型。

BRDF 物理模型是基于电磁波和地表作用后的实际物理过程建立的二向反射模型，物理模型输入的参数都具有明确的物理意义，其优势在于能够模拟冠层内部的多次散射过程。根据物理过程描述方式的不同，又可以将其细分为辐射传输模型、几何光学模型以及二者结合的混合模型。辐射传输模型描述了光子在介质内传播和衰减的规律，代表性模型为Suits 模型（Suits，1971）和SAIL模型（Verhoef，1984）。几何光学模型属于景合成模型，Li和Strahler（1985）创建的Li-Strahler几何光学模型充分考虑了冠层的三维结构对光线的阻拦作用，并由此提出了光照植被、阴影植被、光照背景、阴影背景4 个分量，以这4 个分量的面积比例为权重对四分量的反射率加权求和计算冠层反射率，代表性的模型为GOMS 模型（Li 和Strahler，1992）和四尺度几何光学模型（Chen 和Leblanc，1997）。几何光学模型虽然可以描述树冠的三维结构，但是却无法像辐射传输方程那样刻画树冠内部的多次散射作用，而辐射传输模型正相反。因此，为了充分发挥这两类模型的优势，发展了混合模型，代表性的模型包括GORT 模型（Li 等，1995；Ni 等，1999）、GeoSAIL 模型等（Huemmrich，2001）。近年来，为了满足复杂地表的遥感应用，发展了一批耦合地形的BRDF物理模型。Wen 等（2018）将其分类为单一坡模型（如GOSAILT 模型（Wu 等，2019a）、PLC 模型（Yin等，2017）和GOST 模型（Fan 等，2014）等）和复合坡模型（如dESM模型（Hao等，2018））。

早期，由于对电磁波与地表作用的物理机制认识的不是很清晰，主要靠经验公式对地表二向反射进行模拟。代表性经验模型有Minnaer 模型（Minnaert，1941）、Walthall 模型（Walthall 等，1985）、Staylor-Suttles 模型（Staylor 和Suttles，1986）等。尽管经验模型形式简洁，易于计算，但是依赖大量的实测数据。为了克服经验模型的缺陷，学者们提出了半经验模型来模拟地物的BRDF。半经验模型大多是对物理模型中重要的因子进行参数化而得到的，不仅具有一定的物理意义，还有简洁易算的优点，代表性模型为核驱动模型（Roujean 等，1992）和PRV 模型（Rahman等，1993）。同样，为了使核驱动模型适用于山地复杂地表，发展了一批耦合地形的核驱动模型。如在GoSAILT 模型基础上发展的KDST 模型（Wu等，2019b）、在dESM 模型基础上发展的TCKD 模型（Hao等，2020）以及在山地辐射传输模型基础上发展的LKB_T模型（Yan等，2022）等。

计算机模拟模型主要是通过计算机模拟电磁波在地物表面及内部的辐射传输过程，从而实现地表二向反射特性的模拟。依据模型算法原理的不同，可分为基于辐射度原理方法和基于光线追踪原理方法。基于辐射度原理方法的代表性模型有DIANA 模型（Goel 等，1991）、RGM 模型（Qin和Gerstl，2000）、Rapid 模型（Huang 等，2013）等；基于光线追踪原理方法的代表性模型有DART模型（Gastellu-Etchegorry 等，1996）、FLIGHT模型（North，1996）、Raytran 模型（Govaerts 和Verstraete，1998）、Less模型（Qi等，2019）等。

3 BRDF反演问题

3.1 病态反演问题

Tarantola（1987）指出地学中的参数反演问题，往往都是欠定（Underdetermined）和病态的，即当观测数目不足或者观测数据之间的相关性较高时，参数的解无法固定。同样，在地表BRDF反演过程中也面临着病态反演问题（李小文 等，1997）。

BRDF 表征的是地表在不同太阳—观测角度组合下的反射特征，所以要想反演地表BRDF就需要不同太阳—观测角度组合下的反射率数据。目前应用较广的卫星传感器，如多角度观测的POLDER、MISR，以及宽视场观测的MODIS 都具有多角度观测地表的能力，但是由于扫描方式和轨道高度的差异，传感器获取多角度观测数据的能力和方式也有所不同，如POLDER 相较于MODIS 有着更大的视场角，所以POLDER 能观测到更大天顶角处的反射特征；另外，POLDER 可以从多达15 个角度同时观测地表，而MODIS 等拟多角度卫星只能通过多天观测数据积累才能获取多角度数据集，对于地表特征变化较快的区域，基于此类卫星数据反演BRDF将不再合适。尽管已有众多传感器具备多角度观测能力，但是对于常年多云多雾的山地区域，观测数据不足的问题依然突出，当观测数目小于待反演模型参数时，会导致观测数据的信息量不足，进而造成参数解无法固定的病态反演问题。

除了观测数目不足会导致病态反演问题，观测数据相关性较高也会引起该问题。在BRDF反演中体现为角度分布范围较小，较为集中。首先，大多数的多角度传感器只能在固定或较小太阳入射角变化范围内获取多角度观测数据，如对POLDER BRDF 数据集统计发现，在395 个数据集中有56 个数据集太阳角度采样范围小于5°，约占14%（李小文 等，2000）。另外，由于大气噪声的影响，如云层遮挡可能导致有效观测数据的观测角度也分布在较小范围内。无论是太阳入射角还是观测角度分布集中都会造成观测数据高度相关，数据高度相关意味着多个角度观测数据的信息量和其中一个角度观测数据的信息量相当，但是引入了多个角度观测数据的噪声，这种情况下反演结果受观测噪声的影响很大，精度很低。

3.2 先验知识引入问题

李小文等（1997）认为解决病态反演问题的最好方法就是引入先验知识。学者们也意识到了先验知识的重要性，如Strugnell 和Lucht（2001）收集的68 个地表BRDF 数据集广泛的应用于了BRDF 反演中。但是，需要指出的是引入的先验知识在增加数据信息量的同时也带来了误差和不确定性。先验知识对反演结果造成的误差和不确定性源于多方面原因。（1）先验知识是通过观测统计收集而来的，所以先验知识自身就携带观测噪声，即使先验知识是“完美”的，不含噪声，由于相同地物的二向反射特征在不同时间、空间维度上有差异，仍然无法避免将误差引入到反演结果中。（2）有些BRDF 模型反演约束的先验知识（如核驱动模型的核系数）无法直接通过观测获取，积累收集的先验知识需要利用直接观测的反射率数据转换后才能作为先验知识引入，该过程中会引入转换模型的误差。（3）很多学者在运用先验知识时，通常假设待反演参数的先验分布式服从正态分布，然而研究发现这种假设并不符合实际情况（陈霞 等，2007），所以模型参数正态分布假设同样会给反演结果带来误差。

由于引入先验知识的同时会引入误差和不确定度，所以在引入先验知识之前，有必要权衡引入先验知识后增加的信息量和增加的不确定度，以便确定是否引入先验知识（Wen 等，2017）。在确定引入先验知识后，Li 等（2001a）指出还需要确定引入的先验知识在反演数据中占的比重，因为该参数决定了反演的成败，比重越小，先验知识对方程解的约束作用也就越弱，反之越强。

3.3 观测噪声问题

实际反演过程中，即使观测数目充足，如果忽视观测噪声的影响，也有可能造成反演失败，因为观测数据中可能存在异常值。对于鲁棒性较差的算法，这些异常值可能左右反演的结果，所以反演的第一步就是要剔除观测异常值。在反演过程中，对于越过了“硬边界”的异常值很容易剔除，如反射率值不在（0，1）范围内的就是异常值，可是对于接近“硬边界”的观测数据，无法界定其是否为异常值，如热点方向观测的反射率值会明显高于其它方向的观测值，如果直接剔除可能造成观测数据信息量的损失，这种情况下，我们需要定量刻画不同方向观测数据的噪声水平，即依据噪声水平给观测数据赋权重。

在实际产品生产中，有时为了达到提高运算效率的目的，常假设观测噪声服从标准正态分布。例如，应用于多种BRDF反演算法的普通最小二乘法OLS 就假设观测噪声是服从均值为0，协方差矩阵是单位阵的标准正态分布（Strahler 等，1999；Lucht 等，2000），该假设默认不同观测之间相互独立，且观测噪声水平相同，这样的假设使得在反演过程中不同观测条件下观测数据的权重相同。但是现实情况中，大气状况变化较快，所以每次观测所受的影响也是不同的，晴朗天空下观测数据噪声水平会明显低于多云多雾天空下观测数据噪声水平，所以在反演过程中应该赋予晴朗天空观测数据更大的权重。因此假设观测噪声服从标准正态分布有悖于常识，会给反演结果带来较大误差（Zeng等，2016）。

4 BRDF反演方法概述

反演BRDF的关键在于利用观测数据拟合模型的参数，即方程求解问题，需关注方程是否有解且解是否稳定。当观测数据充足时，可通过最小二乘方法（LSE）求解方程，如对于BRDF 核驱动模型，有研究表明，当观测数目大于7且角度分布范围较大时，可以通过LSE 方法获得稳定解（Lucht 等，2000；Privette 等，1997）。但是，LSE方法假设了观测噪声服从正态分布，而且其鲁棒性较差，受异常观测值的扰动较大。当观测噪声不服从正态分布且观测数据中存在异常值时，可以利用稳健估计方法来代替LSE 方法（Susaki 等，2004）。

当观测数据不足时，或者观测角度分布较为集中时，由于观测信息不足，LSE方法和稳健估计方法反演结果受观测噪声的影响较大，甚至会导致反演失败。针对该问题，学者们提出了不同的解决方案。李小文等（1997）和高峰等（1998）基于模型参数不确定性和敏感性矩阵提出了多阶段目标决策方法，反演过程中通过控制信息流向以充分利用有限的观测数据反演关注的参数。高峰等（2001）提出了“最小方差方法”。不同于最小二乘方法以最小残差和为代价函数，“最小方差法”以白空反照率的最小方差为代价函数。Li 等（1998）认为通过引入先验知识可有效改善病态反演问题，并将先验知识分为3 类：（1）（宽边界形式）待反演参数的取值范围；（2）（δ边界形式）模型中某些参数的值；（3）（软边界形式）参数的联合概率密度分布。随后Li 等（2001a）通过贝叶斯估计方法将软边界形式的先验知识引入到了BRDF 模型反演中，显著提高了反演精度。除了贝叶斯方法，Bell（1978）提出的正则化方法也可以用来改善BRDF 病态反演问题。已发展的BRDF 正则化反演方法可分为正则化平滑法和截断奇异值分解方法。正则化平滑法通过构建平滑矩阵约束最小二乘解（Quaife 和Lewis，2010；Zobitz 等，2020），平滑矩阵的构建也可以视为先验知识的引入。截断奇异值分解方法通过对核矩阵做奇异值分解然后截断极小奇异值获得稳定的方程解，过程中需要通过先验知识获得正则化参数α，所以该先验知识可归为上述的δ边界形式。上述两类方法利用数学技巧改善了BRDF病态反演问题，除此之外，还可以通过扩展观测数据源的手段达到该目的。传统的BRDF 产品多是基于单一的传感器数据生产的，单一传感器对于地面像元的观测次数及观测角度都是极其有限的，为了增加观测次数及扩展观测角度分布，可以联合多源传感器反演BRDF（Samain 等，2006；Liu 等，2010；Wen等，2017）。

然而在某些极端情况之下，可能得不到观测数据。为了应对这种极端情况，学者们将Kalman滤波方法应用到了BRDF 反演中（Qin 等，2006；Samain 等，2008）。对于没有观测数据的窗口，直接利用预测模型对模型参数进行可靠预测；对于有观测数据的窗口，利用观测数据修订预测结果以提高反演精度。

上述反演方法可概括为3 类：基础反演方法、正则化约束反演方法、信息量控制及信息量扩增反演方法，表1总结比较了这些方法的优缺点。基础反演方法无需引入先验知识，最小二乘法和稳健估计方法仅适用于观测数据充足的情况，而最小方差法还适用于观测数据不足的情况。正则化约束反演方法、信息量控制及信息量扩增反演方法需要引入先验知识，主要用于解决BRDF病态反演问题。其中多阶段目标决策方法通过控制观测数据信息流向，让有限的观测信息尽量分配到目标参数反演中。贝叶斯估计法、Kalman 滤波法和多传感器联合反演方法通过扩增数据源的方式解决观测信息量不足的问题。由于方法的原理部分涉及公式推导，所以下面以核驱动模型为例对各方法做详细介绍。

表1 BRDF反演方法特点与适用性总结Table 1 Characteristics and applicability of different BRDF inversion methods

5 BRDF反演方法

5.1 基础反演方法

5.1.1 最小二乘法

式（2）为线性核驱动模型，r表示二向反射率因子BRF；kvol(ti，tv，φ)、kgeo(ti，tv，φ)表示体散射核和几何光学核；fiso、fvol、fgeo表示均一核系数、体散射核系数和几何光学核系数。λ表示波长，ti、tv表示太阳天顶角和观测天顶角，φ表示相对方位角。如果有n组观测数据，核驱动模型可表示为矩阵形式。

式中，M是n维的向量，表示n个不同观测几何下的反射率；K表示核矩阵；X表示核系数矩阵。

最小二乘方法（LSE）方法基于观测数据和模拟数据之间的残差平方和构建代价函数。

式中，Cn表示观测噪声协方差矩阵。对代价函数求最小值，可得核系数的解析解。

LSE是经典的数学回归方法，在地学参数反演中应用极为普遍，但是LSE 方法仅适用于观测数目大于待反演参数的情况。在运用LSE 方法时，尤为需要关注观测噪声的影响，因为LSE 方法的鲁棒性较差，拟合结果易受观测噪声扰动。所以在实际反演中，首先需要剔除异常观测值，然后在此基础上构建观测值权重矩阵，即上式中的Cn。目前提出的构建观测值权重矩阵的方法多是基于像元的二向NDVI 值发展而来的（Gao 等，2002；Li 等，2001b；Zeng 等，2016）。该类方法可以较为准确的为不同观测值分配权重，但是其适用情况有限，仅适用于植被覆盖地表，不适用于非植被覆盖地表，如水体、冰雪覆盖地表以及城市地区。

5.1.2 最小方差法

不同于最小二乘方法以最小残差和为代价函数，“最小方差法”以白空反照率的最小方差构建代价函数。白空反照率方差大小主要取决于观测噪声和BRDF模型自身误差，所以最小方差法拟合了对给定观测值具有最小噪声敏感性的模型参数。核系数X的协方差矩阵可表示如下：

对BRDF在观测方向和入射方向双半球积分后可以得到白空反照率Wb，，Hi是对核函数在双半球积分所得。

式中，θi和θv表示太阳天顶角和观测天顶角；ϕ表示相对方位角。则白空反照率的方差可表示如下，W是由Hi组成的向量。

当观测数目充足时，可利用最小二乘方法（LSE）选择核函数，但是当观测数目不足时，LSE方法得到的结果不够稳定，可靠性不足，且算法鲁棒性较差，易受异常值干扰（Wanner 等，1995）。相较之下，最小方差法有着更好的表现，其不仅适用于观测数据不足的情况，而且在异常值的干扰下，能表现出较强的鲁棒性。Gao 等（2001）借助Kimes 等观测收集的地面BRDF 数据集比较了LSE 方法和最小方差法。实验结果表明，当观测数据充足时，两种方法计算得到的WSA 精度相近。但是当观测数据不足或者观测角度集中时，最小方差法反演的结果更加稳定，可靠性更好，明显优于LSE方法。

5.1.3 稳健估计法

稳健估计方法RE（Robust Estimation）是指当粗差不可避免的情况下，选择适当的方法抑制粗差。其主要思路与最小二乘法（LSE）类似，都是给噪声大的观测值赋小的权重。与LSE 方法不同的是，LSE方法假设了观测噪声服从高斯分布，而RE 方法不需要有此假设，所以RE 方法的鲁棒性更好。本文主要介绍M 估计和最小中值平方法（LMedS）两种稳健估计方法。

（1）M估计。M估计方法通过对模型参数求极大似然估计进行参数估计。核系数的极大似然估计可表示如下：

对于不同的噪声分布形式，基于极大似然估计导出的代价函数估计量也不同，如当观测噪声ε独立同正态分布，极大似然估计可导出最小二乘估计量；当噪声ε独立同双指数分布，如式（10）所示，极大似然估计可导出L1范数估计量，式中，σ表示方差。

对于M 估计方法，可构建一个如式（11）的估计量，其中w(u)表示观测权重，和LSE 方法观测权重构建方法相同，参数s是为了归一化观测噪声。很明显LSE 和L1 范数估计分别是Ψ(u)=u2，|u|的特殊情况。因此如果已知观测噪声的分布形式，将很容易构造出一个稳健估计量。但是，在实际情况中很难知道观测噪声的分布形式，常用的方法是通过选取一些典型的Ψ(u)代入求解，然后比较选择出最合适的Ψ(u)形式用于参数反演。M 估计方法有较高的计算效率，但是其崩溃点较低，即当观测数据中异常值较多时，该方法估算精度较差。

（2）最小中值平方法（LMedS）。LMedS 方法以观测噪声平方的中值作为代价函数进行参数估计。

相较于M 估计方法，LMedS 方法有着较高的崩溃点，能达到0.5，这意味着即使观测数据中异常值比例达到50%，LMedS方法的参数估计依然保持可靠的精度。但是该方法收敛较慢，计算复杂度较高，可借助蒙特卡洛方法提高运算速度。

Susaki 等（2004）评价了普通最小二乘方法（OLS）、M 估计和LMedS 这3 种方法拟合核驱动模型核系数的能力，通过对野外实测的BRF 数据添加噪声生成实验数据。其中M 估计方法用到的是Tukey’s biweight 函数（Huber，1981）。实验结果表明M-estimator 方法和OLS 方法估算的结果相近，LMedS 方法有着更高的估算精度和更好的鲁棒性，其估算的核系数受异常数据的影响最小。

5.2 正则化约束反演方法

5.2.1 Tikhonov正则化平滑法

当观测数目不足或者观测角度分布较集中时，如果直接用最小二乘方法反演模型参数，得到的结果受观测误差的影响很大。为了减弱噪声的影响，可以以平滑矩阵作为正则项约束方程解。通常平滑矩阵可基于先验知识构建，在数学公式上的体现为在最小二乘代价函数的基础之上添加了正则项。

式（14）为最小二乘法的代价函数，式（15）为添加了正则化项的代价函数，其中D表示平滑矩阵，α是正则化系数。对代价函数求最小值，可得核系数的解析解。

式中，平滑矩阵D和正则化参数α是两个关键参数，平滑矩阵D可理解为对方程解约束的形式，可根据核系数随时间变化特征来构建（Quaife 和Lewis，2010；Zobitz 等，2020）。而正则化参数α决定了平滑矩阵约束作用的强弱，作用和贝叶斯估计方法中的先验信息比类似。α的取值范围在（0，1），既不能太大，也不能太小，太大导致过正则化，求得的解远离真值；太小又会导致欠正则化，达不到约束方程解的作用。常用的正则化参数确定方法包括差分准则法（discrepancy principle）、伪—最优判别法（quasi-optimality criterion）及广义交叉验证法（GCV）等方法（Hansen，1998）。上述几种方法都是完全从数学角度出发通过考虑极值问题来确定正则化参数，而忽视了正则化参数的物理意义。另外一些学者针对遥感反演的特点，从信息分析角度出发基于最大熵理论提出了正则化参数确定的方法（杨华 等，2003；赵红蕊 等，2007）。实际生产中，基于式（15）求解核系数的效率较低，为了提高运算效率，可以通过奇异值分解（SVD）方法简化计算过程。

5.2.2 截断奇异值分解法

截断奇异值分解法以模型参数的l2范数构建代价函数。

具体实现方式是通过对核矩阵K进行奇异值分解，然后截断其中较小的奇异值，从而获得稳定的方程解，详细过程如下。首先对核矩阵K进行奇异值分解。

式中，U和V均为单位正交矩阵，U称为左奇异矩阵，U=［u1u2，…，un］，V称为右奇异矩阵，V=［v1v2v3］，Σ=diag（δi），δi（i=1，…，p）是矩阵K的奇异值。那么核系数X可表示如下：

式中，p表示矩阵K的秩，p=min（n，3）。由于M=KX+ε，上式还可以分解为式（17）。

式中，Xtrue表示待反演核系数的真值，ε表示观测噪声。从式（21）不难发现，当处于分母位置的奇异值σi很小时，观测噪声对反演结果的影响很大。为了消除极小奇异值对反演结果的影响，可以将正则化滤波函数χ(δi，α)引入式（21）中，利用滤波函数的滤波作用保留较大奇异值项而淘汰极小奇异值项，在此基础上反演得到稳定的方程解。

在该模型中，最重要的参数是正则化参数α，因为该值决定了哪些奇异值被舍弃掉。该参数可以通过多种方法确定，与4.2.1 节提及的正则化参数选取方法类似，有时为了简单起见，正则化参数还可以直接取传感器波段信噪比倒数SNR-1（Cui 等，2012）。奇异值分解方法不仅适用于观测数据充足的情况，当观测数据不足时，其依然能表现出较好的反演能力，优于Tikhonov 正则化方法和奇异值分解方法（SVD）（Wang等，2007）。

5.3 信息流控制及信息量扩增反演方法

5.3.1 多阶段目标决策方法

多阶段目标决策方法（USM）由李小文等（1997）提出。准确的说，USM 属于一种反演策略，通过控制观测数据信息流向，让有限的观测信息尽量分配到目标参数反演中。其将多参数反演过程分解成多个阶段，每个阶段通过计算不确定性和敏感性矩阵确定该阶段待反演参数以及对待反演参数敏感的多角度数据，具体到反演上还需要配合最小二乘等方法使用。反演过程中需要保证每个阶段反演的参数数量最多不超过4个，循环迭代计算敏感性矩阵、确定待反演参数及敏感多角度数据、基于最小二乘等方法反演参数这一过程，直到收敛。收敛的条件是反演结果满足应用需求，或者观测数据的信息已被充分利用，继续反演下去得到的结果相较于上一阶段的反演结果没有明显变化。李小文等定义敏感性矩阵如下。

式中，ΔBRDF(i，j)表示在第i个采样方向上，当其他参数都固定为期望值而第j个参数在其不确定范围内变化时，模型模拟得到的最大BRDF值与最小BRDF 值的差值。“不确定范围”建议为期望值左右一个标准差范围内。BRDFexp表示第i个采样方向上，所有参数都取值为期望值时，模型模拟的BRDF值。S(i，j)越大，表明第j个参数对第i个方向越敏感，在反演第j个参数时就必须用到第i个方向的观测数据。另外反演时遵循敏感的参数先反演的原则，即S（i，j）越大，反演的顺序越靠前。

每个阶段可以通过选择不确定度和敏感性（US）大的反演参数和多角度数据来固定参数和分割数据集，达到控制信息量的目的。但是当参数的US 较小时，很难确定参数选择和数据集分割方案，基于Shannon熵减理论可以结合USM和信息矩阵，以定量最优控制反演过程中的信息量，从而确定参数选择和数据集分割的最优方案（杨华 等，2003）。

需要强调的是在应用USM 方法时，对于敏感性较强的参数一次反演的结果往往精度不高，需要多次反演逐步提高反演精度，在这个过程中，可以将前一阶段的反演结果作为先验知识加入到后一阶段反演中，能够显著提高参数的反演精度（赵祥 等，2006）。如果待反演模型是非线性模型，参数反演结果将不服从高斯分布。为了将参数反演结果作为先验知识应用于下一阶段的反演，可用左右两个高斯分布函数来表示参数反演结果。此时，“不确定范围”的低边界是期望值和左边高斯函数标准方差之差，高边界是期望值和右边高斯函数标准方差之和。

5.3.2 贝叶斯估计法

贝叶斯估计方法通过引入参数的先验概率分布计算参数后验概率分布，计算得到的最大后验概率对应的参数值即为最优估计值。

式中，X是待估计的参数；P（X）是参数X的先验概率分布；P（M|X）表示参数X确定时，函数值M的概率分布；分母表示全概率，是一个定值。

对于核驱动模型，假设核系数先验概率分布服从均值为X0，协方差矩阵为Cx的高斯分布，则可得到核系数先验概率分布P（X）。

式中，v表示X的维度，det（Cx）表示矩阵Cx的行列式。

假设观测数据服从均值为KX，协方差矩阵为Cn的高斯分布，那么可得到观测数据M的概率分布P（M|X）。

联合式（24）、式（25）和式（26），可求得核系数的最大后验概率函数L（X|M）。

对于式（27），可构建如式（28）所示的代价函数。

Li 等（2001a）首次将贝叶斯估计方法应用到了BRDF 反演中，而后凭借其广泛的适用性成为MODIS BRDF/albedo 全球数据产品生产的基础算法。贝叶斯方法以核系数的先验概率分布这一“软边界”来约束方程解，不仅能解得核系数最优解，还能得到核系数的后验概率分布，有效的改善了病态反演问题。但是在全球范围的产品生产中，如果按照上述方法构建最大后验函数然后求取最大后验概率对应的参数值将耗费大量时间。为了提高效率，可以将核系数先验知识从参数空间（parameter space）转换到数值空间（data space）。具体做法是先利用先验知识的协方差矩阵产生模拟观测数据的核矩阵Ksimu［3］，然后基于先验知识的均值正向模拟观测数据，最后联合真实观测数据反演核系数。在该过程中，需要控制加入的模拟观测数据的比重，该比重可根据已积累的地表BRDF 数据集多次实验得到（Jackson，1979）。需要强调的是在引入先验知识之前必须定量评估先验知识的误差和信息量，因为引入先验知识的时候同时引入了先验知识携带的误差，如果先验知识的误差值大于信息量值，那么对反演结果只会有害无益。

5.3.3 Kalman滤波方法

Kalman 滤波方法属于数据同化方法，对于缺失观测数据的窗口，可以直接通过前一窗口的反演结果预测该窗口参数值；对于有观测数据的窗口，可以利用观测数据对预测值进行修正。具体步骤如下。

（1）预测n时刻待反演参数。基于n-1 时刻的反演结果通过预测模型预测n时刻待反演参数值。

（2）计算n时刻预测结果的协方差矩阵。首先根据上述预测模型计算n时刻预测结果的协方差矩阵，该值可描述预测结果的不确定性。

式中，Q表示模型的噪声方差，刻画了用模型模拟真实物理过程带来的不确定性；Cn-1表示n-1 时刻反演结果的协方差矩阵。

（3）引入n时刻的观测值，基于增益矩阵修订预测值。因为观测值通常并不是待反演参数，如在核驱动模型反演过程中观测数据是反射率，而待反演参数是核系数，所以需要通过观测算子关联预测值与观测值。

式中，H表示线性观测算子，在线性核驱动模型中，H=［1，Kvol，Kgeo］；v表示第n次观测的随机白噪声，其方差为R。然后，就可以基于增益矩阵计算n时刻待反演参数的均值和协方差矩阵。

式中，Gn表示增益矩阵，形式如下：

Kalman 滤波方法与贝叶斯估计方法类似，也需要引入核系数的先验概率分布，因为第一次迭代时的均值X0和协方差矩阵C0需要通过先验知识获得。另外，Kalman 滤波方法也能得到核系数的后验概率估计，与贝叶斯估计方法不同的是，当窗口缺失观测数据时，贝叶斯估计方法只能基于先验知识反演核系数，而Kalman 滤波方法可通过预测算子F可靠预测核系数，F多是对参数随时间变化的统计建模（Xiao 等，2011；Samain 等，2008），反映了参数随时间的变化规律。所以Kalman 滤波方法可认为是在贝叶斯估计方法的基础上增加了参数在时间维度上的变化信息，增多了BRDF反演的信息量。正是基于该特点，Kalman滤波方法可用于BRDFalbedo 产品的时空填补，其填补的产品能表现出较好的时间、空间连续性（Samain 等，2008）。Kalman 滤波方法仅适用于系统为线性且观测噪声为高斯白噪声的情况，对于非线性系统，可以考虑集合Kalman滤波方法（Qin等，2006）。

5.3.4 多传感器联合反演

为了增加观测数目和增大角度分布范围，可以联合不同类型的传感器反演BRDFalbedo。在融合不同传感器时，需要考虑光谱融合、方向融合、时间融合和空间融合4个维度的融合。方向融合是指将不同太阳—观测几何下的观测数据融合在一起，核驱动模型正好可以实现这一目的。时间融合是为了将某一时间窗口内不同时间观测数据融合在一起，如MODIS传感器BRDFalbedo产品就是将16 天观测窗内的观测数据融合在一起，时间融合的主要策略是通过对不同时间的观测数据赋予不同观测权重实现的。空间融合的目的是将不同传感器不同空间分辨率的观测数据融合在一起，常用的空间融合方法有基于频率分析的傅里叶变换（Schowengerdt，1980）和小波变化方法（Garguet-Duport 等，1996），这些方法要求影像中不能有空白像元。另外，对于不规则且分辨率不同的传感器数据，通常借助于多尺度Kalman 平滑方法（Chou 等，1994）实现多源数据空间融合。传感器光谱主要由中心光谱波段值和波谱响应函数表征，对于不同传感器这两者是不同的，所以在实现传感器融合时，还需要考虑光谱融合，这里做详细介绍。

光谱融合其实就是要实现光谱的归一化，其核心思路是将波长参数从二向反射率中分离出来。基于该思路，学者们提出了多元线性回归法、多角度和多波段（ASK）模型方法、多传感器联合反演BRDF 方法（MCBI）等方法（Samain 等，2006；Liu 等，2010；Wen 等，2017）。多元线性回归法是以某一传感器作为参考传感器，建立该传感器某一波段的反射率和另一传感器多个波段反射率的线性回归关系，拟合得到转换系数和截距即可以实现向参考传感器波段的转换（Samain 等，2006）。多角度和多波段（ASK）模型方法通过将波长引入几何光学核函数和体散射核函数中，从而实现多源传感器波段融合（Liu 等，2010）。多传感器联合反演BRDF 方法（MCBI）以MODIS 的7 个波段作为基础波段，通过将AVHRR、VIIRS、MERSI 等传感器波段反射率分解为MODIS 的7 个波段反射率以实现波长参数的分离，该方法的关键是通过计算“净信息指数”决定是否融合新的遥感数据（Wen等，2017）。

式中，K表示核函数矩阵，对K′K做Schur 分解，可得特征值对角阵V=diag（λi），λi（i=1，…，21）和特征向量G，然后计算“净信息指数”。

式中，NII 表示“净信息指数”，I表示信息指数，MSE 表示均方误差。若NII 大于0，表示引入数据的信息量大于携带的噪声，允许引入，反之则拒绝引入。

6 结语

地表二向性反射的研究对于地表参数的反演有着重要的意义，经过20 世纪80 年代以来的发展，地表BRDF的反演取得了显著的进展。本文总结了BRDF反演的各种方法，这些方法充分利用发掘了先验知识信息及观测数据信息，缓解了反演过程中因观测数据不足或者观测噪声导致的病态反演问题。在反演BRDF过程中，重点关注了观测噪声及观测异常值、观测数量不足或观测角度分布集中两方面问题。

（1）利用卫星观测数据反演地表BRDF时，通常需要多角度数据，多角度传感器或者宽视场传感器多天积累的观测数据可提供多角度数据。由于受大气的影响，多角度数据的大气校正结果常伴有较大的噪声。对于异常数据，可利用非对称高斯拟合函数（Jonsson 和Eklundh，2002）和分段logistic 拟合函数（Zhang 等，2003）等时间平滑法初步剔除。对于有着不同观测噪声水平的观测数据，在参与反演模型参数过程中需要赋不同权重，在数学公式中的体现为构建观测噪声协方差矩阵。目前常用的方法是基于地面像元的NDVI 通过迭代计算观测数据权重（Gao 等，2002；Li 等，2001b；Zeng 等，2016）。当观测噪声不服从高斯分布时，还可以借助于稳健估计方法估算模型参数。

（2）观测数据不足或观测角度分布集中时，会导致观测信息量不足，进而影响BRDF 反演精度。为了减小因信息量不足而导致的病态反演问题的影响，本文总结了多种方法，如最小方差法以最小噪声敏感性构建代价函数，Tikhonov正则化平滑法通过平滑矩阵约束方程解，截断奇异值分解法通过剔除极小奇异值来抑制误差影响等。多传感器联合反演方法通过扩展数据源的方式增加信息量，贝叶斯估计方法、Kalman 滤波方法通过引入参数的先验概率分布来增加信息量。多阶段目标决策方法基于不确定性和敏感性矩阵控制观测信息量流向，以发掘有限的观测信息反演关注的参数。

虽然目前针对地表BRDF遥感反演方法和策略的研究取得了较好成果，然而要想反演得到高质量的BRDF产品，最重要的还是获取高质量、充足数量的观测数据以及适合地表BRDF反演的遥感模型。随着国内外不同尺度遥感卫星数据的丰富，对陆表BRDF遥感反演也提出了更高的需求，但仍存在着如下问题：

（1）多尺度多地类BRDF 模型体系仍未建立，从现有可支持遥感反演的BRDF模型看，主要适用于中低分辨率遥感像元尺度，普遍缺乏对高分辨率像元尺度的支持，这是由于高分辨率像元间具有较强的邻近效应及相互遮挡效应，由中低分辨率像元发展而来的BRDF模型不再适用于高分辨率遥感数据，从而导致高分辨率BRDF反演困难。随着国内外高分辨率卫星技术和无人机技术的快速发展，发展适合高分辨率像元尺度的BRDF模型及其反演方法已迫在眉睫。

（2）支持BRDF遥感反演的多角度数据普遍不足，BRDF 先验知识缺乏。BRDF 的高精度遥感反演始终离不开多角度数据的支持，然而由于天气等因素的影响，多角度数据普遍缺乏，导致当前的BRDF反演过程中合成周期较长，特别是在山区复杂天气情况下，无论是中低分辨率数据还是高分辨率数据，其有效观测的多角度数据更加缺乏。即使引入BRDF 先验知识可提高BRDF 反演的精度，但现阶段山区地表和高分辨率BRDF先验知识同样缺乏。因此，发展多源、多尺度数据联合反演方法以及构建山区地表和高分辨率地表BRDF先验知识是实现均质复杂地表和多尺度像元地表BRDF反演精度提高的关键所在。

（3）过去几十年，随着遥感数据的积累，遥感已进入了“大数据时代”。遥感大数据的分析应用对于全球范围内自然灾害的监测、精准农业的实施等有重要意义。目前全球规模最大的地球观测组织GEO（Group on Earth Observations）已收集了来自全球超过80 个国家的7000 万条观测数据，各国也在建立自己的遥感大数据体系。面对海量的“遥感大数据”，人工智能技术成为了处理分析“遥感大数据”的重要手段。近年来，人工智能技术在遥感影像变化监测（Shi 等，2020）、遥感信息挖掘（Wang 等，2020）等方面得到了较好的应用，并成功应用到了大范围区域土壤水分监测（Sahoo 等，2017）、粮食估产（Wolanin 等，2020）等领域。“遥感大数据”同样为BRDF 反演提供了一个巨大的先验知识库，但是如何利用人工智能技术从规模巨大的多源异构“遥感大数据”中发掘并构建BRDF先验知识库值得进一步研究。