地表反照率（Albedo）定义为地面反射的太阳辐射与到达地面的太阳辐射之比，它代表了晴空无云地球表面对太阳辐射的反射能力, 是数值模拟和地面能量平衡方程中的一个重要参数^[1]。地表反照率的确定对研究下垫面的热力性质、地表能量平衡和全球气候变化有重要意义。准确地测定地表反照率是地表辐射平衡和气候研究中的一项重要工作^[2]。

由于反照率是反射率在入射和出射半球的积分，因此需考虑地表反射率的各向异性特征。目前的研究中主要采用双向反射分布函数（bidirectional reflectance distribution function，BRDF）描述地表反射率的各向异性特性。BRDF为入射方向地表辐射微增量与其所引起的反射方向反射辐亮度增量间的比值^[3-4]。为了使BRDF模型更具有真实性，多数BRDF模型介于镜面反射与漫反射之间，并设定BRDF具有3个基本性质：微元性、能量守恒性、霍茨互异性^[5-6]。现在基于遥感数据反演反照率的模型主要包括物理模型（辐射传输模型、几何光学模型、混合模型）、经验模型（Minneart、Walthall模型）、半经验模型（核驱动模型、Hapke模型）、计算机模拟模型（蒙特卡罗模型）等^{[4, 8-10]}。此次反演BRDF将会用到半经验核函数模型^[8-11]。最经典的核驱动BRDF算法已经十分成熟^{[4, 8-11]}，即先进行云检测与云去除、大气校正、几何矫正，然后挑选BRDF核驱动模型进行回归，最后进行BRDF半球积分和窄波段向宽波段转换^{[9, 12]}。

尽管如此，反照率反演依然面临一系列问题，如对数据量少、角度数据不全的卫星进行反演时可能无法反演或精度较低；对时间分辨率低、时序数据不全情况进行处理，时间过长可能面临植被生长、融雪等问题；再如地形的变化对反照率的影响等。近年来遥感科学家提出了多种处理方法，如用正则化滤波约束^[13-15]或先验知识算法（如角度格网法（angular bin，AB1））反演数据不全的卫星^{[2, 9]}，通过多卫星多传感器相互补充综合反演解决时序过长的问题等；并且建立了相关模型去矫正地形的影响，如邻坡反射影响等^[16]。这些方法中，正则化滤波方法能较好地优化角度数据不足及观测样本不充分下BRDF反演欠定问题^[13-15]。而在观测角度单一的情况下，AB1算法能很好地弥补BRDF算法的缺陷。

高分四号卫星（GF-4）为地球同步轨道卫星，缺乏多角度观测数据。因此本文用AB1算法进行GF-4地表反照率反演，以求进一步优化和验证AB1算法的准确性，对算法做出一些改进与讨论；同时针对GF-4观测角度数据的不足，改进得到高质量的产品，将有助于我国带动行业部门推广应用高分数据，推动高分应用业务的发展。试验结果表明，AB1算法十分适合GF-4数据，能得到较好的反照率结果。

1 反演方法

GF-4在高度3.6万公里的地球同步轨道运行，采用面阵凝视方式成像，可观测面积大，还能长期对某一地区固定观测，各波段详细信息见表 1。

地球反射极化和方向测量仪（POLarization and Directionality of the Earth’s Reflectances，POLDER）传感器是由法国空间中心研制并搭载在法国小卫星PARASOL（Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a LiDAR）上的传感器，是世界上第一个可进行多角度偏振成像观测的卫星载荷。POLDER传感器包括从蓝波段(0.443 μm)到近红外波段(1.020 μm)的9个通道，分辨率为15~30 m，星下分辨率约为5.3×6.2 km^[12]，偏振观测在0.490 μm、0.670 μm和0.865 μm，光谱波段的带宽在20~40 nm之间^[17]。

1.1 BRDF反演过程

BRDF函数即反射率关于太阳天顶角、太阳方位角、观测天顶角、观测方位角4个角度的函数，基于BRDF关于主平面对称的假设(互异性)，太阳天顶角和太阳方位角通常可由一个相对方位角代替（即二者差值的绝对值）。由于AB1算法中对POLDER数据进行先验处理时需要利用BRDF反演得到反照率，因此可将BRDF反演视为AB1算法的基础。BRDF反演过程主要分为BRDF拟合、BRDF球面积分、窄波段转宽波段3步。

1) BRDF拟合

本文采用基于半经验核驱动BRDF模型来估算窄波段反照率，选择核驱动BRDF函数如下^[18-22]：

$ \begin{array}{l} \begin{array}{*{20}{l}} {R\left( {\theta ,\vartheta ,\varphi } \right) = {f_{{\rm{iso}}}} + {f_{{\rm{vol}}}} \cdot {k_{{\rm{vol}}}}(\theta ,\vartheta ,\varphi ) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {f_{{\rm{geo}}}} \cdot {k_{{\rm{geo}}}}(\theta ,\vartheta ,\varphi )} \end{array}\\ \end{array} $

(1)

式中，R(θ, ϑ, φ)是某一像元的单波段反射率，是观测几何的函数；θ、ϑ、φ分别为太阳天顶角、观测天顶角、相对方位角^[18-19]；f_iso、f_vol、f_geo是待定系数，分别表示各向同性核系数、体散射核系数和几何光学核系数^[22-23]，即需要拟合出的目标；k_vol、k_geo分别是体散射核函数和几何光学核函数^[20-22]。本文选择体散射核为Li-SparseR^[4]、几何光学核为Ross-Hostpot的核函数^[18-22]。

因此，对于一个像元，只要有3组以上的观测数据（即反射率R以及对应的观测几何），就可以按照式（1）用最小二乘法拟合出该像元的BRDF函数。

2) BRDF球面积分

自然状态下，下行辐射通量可分为直射光和散射光两个部分，地面对这两部分的反射特性有所不同。本研究拟用黑空反照率(black-sky albedo，bs)和白空反照率(white-sky albedo, ws)来描述地表反照率特性。黑空反照率又称方向-半球反照率，即将BRDF在反射半球进行积分，表示仅有太阳直接入射的地表反照率，是太阳角度的函数：

$ {\alpha _{{\rm{bs}}}}(\theta ) = \smallint _{0{\rm{ }}}^{2\pi }\smallint _0^{\frac{\pi }{2}}R(\theta ,\vartheta ,\varphi ){\rm{sin}}\theta {\rm{cos}}\theta {\rm{d}}\theta {\rm{d}}\varphi $

(2)

白空反照率又称漫射半球-半球反照率，即将BRDF既在反射半球积分，又在入射半球积分，表示仅有大气漫散射的地表反照率：

$ {\alpha _{{\rm{ws}}}} = \smallint _{0{\rm{ }}}^{2\pi }\smallint _0^{\frac{\pi }{2}}{\alpha _{{\rm{bs}}}}(\theta ){\rm{sin}}\theta {\rm{cos}}\theta {\rm{d}}\theta {\rm{d}}\varphi $

(3)

最终，波段λ的地表蓝空反照率近似等于白空反照率和黑空反照率的加权组合：

$ \alpha (\theta ) = {\rm{ }}(1 - S){\alpha _{{\rm{bs}}}}(\theta ) + S{\alpha _{{\rm{ws}}}} $

(4)

式中，S是权重系数，取决于大气状况，利用大气辐射传输模型（second simulation of the satellite signal in the solar spectrum，6S）模型，建立S与气溶胶光学厚度的函数关系，从而确定S的值，其中气溶胶光学厚度来自MOD04产品。因此，将BRDF函数进行半球积分后，就得到窄波段白空、黑空反照率，加权后得到最终蓝空反照率。

3) 窄波段转宽波段反照率

前文计算得到了单个波段的反照率，但最终需要得到宽波段的反照率。宽波段反照率是反照率关于波长的积分。在没有函数关系，只有几个窄波段反照率（绿、蓝、红、近红外）的前提下，对于每一类地物（分类标准如植被、水体、雪地、裸地等），宽波段反照率与多个波段的窄波段反照率之间存在线性经验关系^[2]:

$ \begin{array}{*{20}{l}} {\begin{array}{*{20}{l}} {F(x) = {a_0} + {a_1} \cdot f({x_1}) + {a_1} \cdot f({x_1}) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {a_2} \cdot f({x_2}) + \cdots + {a_n} \cdot f({x_n})} \end{array}} \end{array} $

(5)

式中，F(x)为宽波段反照率；f(x₁)，f(x₂)…f(x_n)为第1，2…_n个波段的窄波段反照率（如红、绿、近红外等）；a₀，a₁，a₂…a_n为波段转换对应的待求系数。

为了求得适用于POLDER传感器的经验关系系数，需要利用POLDER传感器相应_n个窄波段通道的光谱响应函数，以及ENVI(the environment for visualizing images)所附带的190条典型地物光谱曲线。光谱响应函数描述了窄波段范围里传感器对每一个波长值的响应率。典型地物光谱曲线则描述了每一种典型地物在每一个波长值的反射率。将光谱响应函数与相应窄波段范围内的光谱曲线做卷积，就得到了每个窄波段的反照率。而单独将光谱曲线做加权平均，可以得到宽波段反照率。

将190条曲线所得到的宽波段反照率与卷积得到的窄波段反照率进行回归，就可以得到非冰雪类（式（6））与冰雪类(式（7）)的经验系数：

$ \ \ A1{\rm{ }} = {\rm{ }}0.283{\rm{ }}27{\rm{ }} \times {a_{{\rm{blue}}}} + {\rm{ }}0.154{\rm{ }}21{\rm{ }} \times {a_{{\rm{green}}}} +\\ 0.312{\rm{ }}95{\rm{ }} \times {a_{{\rm{red}}}} + {\rm{ }}0.250{\rm{ }}67{\rm{ }} \times {a_{{\rm{nir}}}} - {\rm{ }}0.000{\rm{ }}93 $

(6)

$ \ \ A2{\rm{ }} = {\rm{ }}0.556{\rm{ }}77{\rm{ }} \times {a_{{\rm{blue}}}} + {\rm{ }}0.045{\rm{ }}37{\rm{ }} \times {a_{{\rm{green}}}} -\\ 0.220{\rm{ }}79{\rm{ }} \times {a_{{\rm{red}}}} + {\rm{ }}0.354{\rm{ }}94{\rm{ }} \times {a_{{\rm{nir}}}} + {\rm{ }}0.264{\rm{ }}59 $

(7)

1.2 AB1算法过程

AB1算法的核心思想是：宽波段反照率与某一个特定角度的几个波段反射率之间存在一种线性经验关系^{[2, 24-25]}，即:

$ {\alpha _{{\rm{bs}}}}(\theta ) = {c_0} + {c_1}{R_1} + {c_2}{R_2} + {\rm{ }} \cdots {\rm{ }} + {c_n}{R_n} $

(8)

$ {\alpha _{{\rm{ws}}}} = {d_0} + {d_1}{R_1} + {d_2}{R_2} + {\rm{ }} \cdots {\rm{ }} + {d_n}{R_n} $

(9)

式中，R_n是特定角度下第n个波段的反射率；c_i与d_i(i=0, 1, 2…n)是回归系数, 即查找表。如果求出查找表，就算GF-4的角度数据不全，只要给出一个角度下的一景图像，也能求出宽波段反照率。而查找表则需要先通过POLDER卫星的数据集进行计算，然后利用POLDER卫星与GF-4的光谱响应函数之间的差异进行一定的转换，最终得到适合于GF-4的查找表。AB1算法的主要过程也分为3步，即POLDER数据集计算、格网查找表建立、GF-4卫星宽波段反照率计算。

1) POLDER数据集计算

现有的POLDER数据集每一类地物都给出了若干像元不同角度的观测数据。对于每一个像元均可进行BRDF反演，得到BRDF函数以及相应的宽波段反照率。最终对于每一类地物可得到若干像元的BRDF函数与宽波段反照率，以便于求经验关系，构建查找表。

2) 格网查找表建立

首先建立GF-4卫星格网回归查找表，也就是建立GF-4卫星宽波段反照率（白空、黑空）与每一角度下（即一格网）4个波段反射率之间的回归关系:

$ {\alpha _{{\rm{bs}}}}(\theta ) = {c_0} + {c_1}{R_1} + {c_2}{R_2} + {\rm{ }}{c_3}{R_3}{\rm{ }} + {c_4}{R_4} $

(10)

$ {\alpha _{{\rm{ws}}}} = {d_0} + {d_1}{R_1} + {d_2}{R_2} + {\rm{ }}{d_3}{R_3}{\rm{ }} + {d_4}{R_{n4}} $

(11)

对于每一个像元，可通过BRDF函数算出每一角度下的反射率。为了将不同角度反射率离散化，需将POLDER的BRDF进行格网划分。将太阳天顶角按每2°间隔，0°~80°范围的观测天顶角按每2°间隔，0°~64°范围的相对方位角按每5°间隔，0°~180°范围划分成41、36、37个格网。对每一个格网计算该角度值下的4个波段反射率。然后利用POLDER与GF-4的光谱响应函数以及190条典型地物的光谱曲线做卷积，将每个格网POLDER 4个波段反射率转换到GF-4的4个波段反射率。得到对应波段回归关系如下：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{gf}}4{\rm{blue}} = 0.99194 \times {\rm{POLDE}}{{\rm{R}}_{{\rm{blue}}}} - 0.00154}\\ {{\rm{gf}}4{\rm{green}} = 0.99769 \times {\rm{POLDE}}{{\rm{R}}_{{\rm{green}}}} - 0.00227}\\ {{\rm{gf}}4{\rm{red}} = 0.99860 \times {\rm{POLDE}}{{\rm{R}}_{{\rm{red}}}} - 0.00272}\\ {{\rm{gf}}4{\rm{nir}} = 0.99987 \times {\rm{POLDE}}{{\rm{R}}_{{\rm{nir}}}} - 0.00181} \end{array}} \right. $

(12)

此处光谱转换与原AB1算法中有所不同，原算法是采用单波段对多波段的多变量回归，而此处采用的是单波段对单波段回归。单波段转单波段精度会比多波段转单波段略低，均方根误差（root mean square error，RMSE）为0.02左右，本文计算中精度已经足够，相关系数r均大于0.99, 且单波段转单波段最大的优势在于能够保证波段内容的纯粹以及波段间的一致性。POLDER与GF-4 4个波段光谱响应曲线如图 1所示。

图 1由从中国资源卫星应用中心下载的GF-4光谱响应率数据与法国空间研究中心下载的POLDER光谱响应率数据所绘制。从图 1中可以看出，两个传感器每个对应光谱响应曲线的峰值十分吻合，只是波宽存在差异。而每一个像元的POLDER宽波段反照率不需要再利用光谱响应函数进行转换。每一类地物都有若干像元的数据可用于回归。现将宽波段反照率与每一个格网的4个窄波段反射率进行回归，可得到每一个格网里每一类地物的宽波段反照率与窄波段反射率之间的回归关系式（10）、（11），其系数随角度格网的改变以及地物类别不同而不同。此系数表即查找表。

3）GF-4反照率计算

代入一景GF-4的反射率图像即可通过查找表反演出宽波段反照率。在代入图像时，还需给出每个像元的角度数据。太阳天顶角、观测天顶角、相对方位角可通过卫星高度、星下点经纬度、拍照地面点经纬度、拍照时间算得，具体几何求算过程见文献[26]。

2 实验 2.1 研究区介绍

为了充分体现算法的适用性，本文选取河北和内蒙古两个区域，分别用中分辨率成像光谱仪（Moderate-resolution Imaging Spectroradiometer，MODIS）和Landsat 8产品进行算法测试。分别选取北京市部分区域以及河北怀来试验场、内蒙古部分地区作为研究区，如图 2、图 3所示。研究区域包括植被、水体、裸地等多种地物类型，地物覆盖类型丰富。数据分别拍摄于2016-04-10 10:54和2016-08-03 10:33。

2.2 反照率结果验证与分析

由于地面验证数据不全，此处选用MODIS与Landsat相应区域的反照率产品来做交叉验证, 交叉验证也是反照率产品适用性验证的一种常用方法^[27]。选取成像日期与GF-4影像数据的拍摄日期最为接近的MODIS MCD43A3反照率产品，该产品提供了全球8 d合成的500 m地表短波反照率。此次对比了验证区域短波波段的反照率。Landsat 8没有直接反照率产品，因而可以选择反射率基础数据，先利用MODIS相应时空范围先验BRDF产品进行窄波段反照率转换^[28]，然后通过Liang^[29]建立的窄波段反照率向宽波段反照率的转换系数求得。

利用GF-4数据反演的地表反照率结果与MODIS地表反照率产品进行对比，如图 4所示。

在图 4研究区域内，先将GF-4反照率影像升尺度到500 m，然后生成400个随机点提取GF-4与MODIS的地表反照率值进行对比，如图 5所示。从图 5中结果可以看出，GF-4反照率与MODIS反照率一对一散点近似分布在直线y=x附近，且RMSE为0.027 3，优于通常对反照率产品的精度要求（RMSE达到0.05）^[2]。因此可以认为，由AB1算法反演GF-4反照率产品具有较高的精度。

将利用GF-4数据反演的地表反照率结果与Landsat 8地表反照率产品进行对比，如图 6所示。在图 6研究区域内随机挑选400个点，对这400个点利用GF-4数据获得的地表反照率值与Landsat 8数据计算得到的地表反照率值进行对比，如图 7所示。

从图 7中结果可以看出，GF-4反照率与Landsat 8反照率一对一散点近似分布在直线y=x上，线性非常好，且均方根误差为0.017 6。由此可见，由AB1算法反演GF-4反照率与Landsat 8反照率产品线性关系较好，精度较高。

3 结语

本文研究可验证AB1算法在GF-4这类高空间分辨率且缺少多角度观测数据的静止轨道卫星上进行反照率反演的可行性，以利于推广高分系列卫星的应用。同时，本文基于实现AB1算法过程中遇到的困难，对AB1算法的一些细节进行了改进探讨，包括改进光谱转换过程，讨论查找表分类以及产品与验证的时间空间尺度效应等。研究结果表明，利用AB1算法做单一观测角度反照率反演结果有较高精度，RMSE分别为0.027 3与0.017 6，证明了AB1算法思路的合理性，同时也实现了针对GF-4角度数据不足来进行反照率反演的问题，对地球同步轨道卫星的参数反演具有一定的参考价值。

本文还有如下问题值得进一步讨论。

1）空间分辨率问题。AB1算法本身没对空间分辨率做出更多研究。算法的查找表系数只针对了地物类别而没有针对像元尺寸。而实际上构建查找表的过程可能与像元尺寸相关。BRDF模型作为构建查找表的基础具有微元性，像元尺度放大会降低BRDF精度，所以训练数据集应具有高空间分辨率。至于回归过程与测试过程，如果像元尺度过大，可能会导致像元类别不纯、几何角度不清晰等问题，进而影响到精度。因此用AB1算法运算不同分辨率的影像（如GF-4和MODIS）并进行对比值得进一步探讨。

2）时间分辨率问题。传统的BRDF反演需要多角度观测数据，通常需要长时序观测。一些研究通过多源遥感数据协同反演来弥补这一不足^[30]，而AB1算法的单一图像反演特性使得AB1算法可以获取某一时刻的反照率。本文虽然做了两组不同时间空间分辨率的验证，却仅仅初步体现了AB1算法对GF-4影像的适用性，由于数据原因，没有进一步做更多的不同时空分辨率的对比实验。可进一步结合BRDF以及多源遥感数据反演，来研究AB1算法时间分辨率的问题。

3）分类问题。地表分类问题是AB1算法的另一个重要因素，包括查找表系数和光谱转换系数的地表分类。且地表分类的精细程度与空间分辨率之间也有交叉影响。本文曾试图降低查找表分类标准，并与MODIS和Landsat影像进行比较，发现在与高分Landsat影像比较时精度略有下降。如何优化分类标准能够既不繁杂又有代表性地将查找表聚散特征分离开，从而在不同分辨率影像下达到较好精度，需要进一步探讨。并且某些地类如冰雪和水体，在构建查找表时添加一些单独的模型（如后向散射核^{[2, 14]}）效果会更好。

此外，基于本文还可以引入AB2算法，根据相关对比方法^[31]加强验证，结合正则化方法等。