1、引　言

地表反照率定义为地表反射的太阳入射辐射的能量在半球空间的积分与入射太阳辐射能量的比值，它是研究地表能量收支平衡以及全球气候变化的重要参数，目前已被广泛应用于众多的气候模型(Dickinson，1983)。由于地表反射的太阳入射辐射是各向异性的，并且随入射和出射方向而变化，因此精确反演地表反照率需要考虑地表反射的各向异性影响。目前，针对在轨运行的卫星遥感数据以及部分航空遥感数据已有大量地表反照率反演研究，并进行了有效的反照率产品的地面验证(Lucht 等，2000；Schaaf 等，2002；Liang 等，2003；董亚冬 等，2014)。其中，基于线性核驱动Ross-Li模型的地表BRDF反照率遥感反演算法是目前地表反照率遥感反演中应用最广泛的方法。已经在百米到公里级的MODIS、MISR的地表反照率产品中得到了很好的应用(齐文栋 等，2014)。

随着精细化监测的不断发展，中高空间分辨率的地表反照率产品也逐渐成为气候模型的主要输入。中高分辨率遥感具有较好的地物细节信息，有助于人们更好地了解地面的真实情况。然而，受空间采样方式的限制，30 m分辨率的Landsat卫星数据通常只有一个近天顶方向的观测，对其进行地表反照率反演通常需要借助先验知识(Li 等，2001)。Shuai等利用30 m分辨率的Landsat地表反射率以及相同时间高质量的500 m空间分辨率的MODIS BRDF产品作为先验知识，反演了30 m分辨率Landsat地表反照率(Shuai 等，2011)。然而，Landsat反照率至今主要采用特定、离散的地表通量塔的时序反照率观测数据进行验证。由于通量塔站点的设置和维护需要花费大量精力，其数量相对有限。因此，目前针对区域的验证仍是我们面临的一个较大的问题。

机载飞行试验可以快速获取区域高空间分辨率(通常为米级)的多角度观测，为实现星载反照率的区域验证提供了契机。机载红外广角双模式成像仪WIDAS(Wide-angle Infrared Dual-mode line/area Array Scanner)是2016年在根河森林区域开展的大型实验中的主要机载传感器之一，它能够获取5个可见光波段的12个观测角度的反射率数据。但是，由于大角度下的观测数据存在一定问题，多角度观测数据只在25°范围内可用，从而影响了WIDAS多角度数据所包含的信息量。BRDF原型反照率反演算法为有效地解决这一问题提供了新途径。验证表明，该算法能够在小观测角度时实现较高精度的反照率反演(Li 等，2001；Jiao 等，2014；Dong 等，2016)。这一方法已成功应用于2008年和2012年黑河农田区域的地表反照率产品生产和验证(Jiao 等，2012；张虎 等，2013；Jiao 等，2014)。基于该算法和机载数据反演的1 m高空间分辨率反照率(Jiao 等, 2015)，在保证产品本身精度的条件下，为验证目前主流的30 m分辨率的Landsat反照率产品提供了新的思路。

本文首次将BRDF原型反演算法应用于森林区域的机载WIDAS观测数据的反照率反演，并基于WIDAS反照率对Landsat反照率进行了初步验证。探索了机载反照率用于中高分辨率的Landsat反照率产品的验证潜力。

2、数据源     (2.1) 根河区域机载和验证点数据

研究区的内蒙古自治区大兴安岭森林生态系统国家野外观测研究站位于内蒙古自治区呼伦贝尔市根河林业局潮查林场境内，地理坐标为50°49′—50°51′N、121°30′—121°31′E，海拔800—1000 m。植被类型为以兴安落叶松为主的北方针叶林，在原始林通量观测场设置森林气候梯度观测系统，装有总辐射、净辐射、光合有效辐射、空气温湿度等13项71个传感器自动观测记录存贮装置。

2016年8月19日，试验在1 km²范围内，利用运五飞机进行了WIDAS(多角度集成平台)和LiCHy(激光雷达、CCD、高光谱多传感器集成平台)的观测，同时利用中国科学院遥感与数字地球研究所研发的无人机平台进行了高光谱、热红外、CCD和PAR的观测。WIDAS传感器为广角镜头成像，经过数据预处理后，可以从相邻图像的重叠区域提取同一地物的多角度信息，WIDAS参数见表1。飞行区域及路线见图1，黄色区域为1 km×1 km重点试验区，即无人机高光谱/红外作业区。红色矩形为600 m飞行高度LiCHy的飞行区域，黑色矩形为1000 m飞行高度LiCHy的飞行区域，蓝色矩形为1500 m飞行高度WIDAS的飞行区域。本文的研究区域具体范围为1896 m×2712 m。本次航空飞行多角度观测试验区域小且观测角度多，单轨飞行路径短且飞行轨迹是米字型交叉，因此能够获取多个平面的观测。

相对于Landsat星载数据，航空遥感更容易在短时间获取小区域高精度的多角度数据，且受大气影响小，因此，大气校正精度相对较高。

如图2，试验区内有一个验证点，通量塔高64 m，反照率表架设距离冠层顶约24—30 m，每隔10 min采样一次(十分钟内的观测均值)，这为反照率的进一步校验提供宝贵的地表通量塔站点数据。但不同于黑河区域的农作物试验区，根河森林试验区内布置临时反照率表困难。因此，本研究没有获取更多地表站点反照率数据用于验证。

本研究采用飞机过境时间为中心，将前后各10 min内的观测数据的平均值作为地面观测的真值。为了消除噪声的影响以及与地面站点的尺度进行匹配，反照率表的观测范围可以通过下式进行估算

$ f = 2H\tan \left( {{\rm{HFOV}}} \right) $

(1)

式中，H为反照率表的高度，HFOV为半视场角。采用的反照率表(CNR4)的半视场角为85°，通过上式计算得到的反照率表的视场范围约为570 m。为了保证验证数据在观测视角一致，将WIDAS反照率结果在地面站点附近约450 m×450 m内的反照率的均值与地表实测数据进行比较。

    (2.2) Landsat 8 OLI/TIRS数据

对于现代的遥感应用，尤其是局部或区域的土地利用和气候研究，对遥感数据的分辨率提出了较高的要求。在许多的卫星对地观测系统中，Landsat系列卫星由于具有40多年的观测历史，分辨率较高(30—60 m)、免费共享政策等优点，成为最为广泛应用的系列卫星之一。最新的Landsat 8卫星自2013年2月11日发射以来，已在轨运行多年，它以每天550幅的速率向地面传回了海量的对地观测数据，并与Landsat 7 ETM+一起，构成了8 d间隔的Landsat重复观测周期，为地球表面的生态环境变化监测提供了极为宝贵的数据(徐涵秋，2015)。

本研究使用了2015年全国Landsat 8经过大气校正和地形辐射校正的地表反射率高级产品数据(Hu 等，2014；Wang 等，2016)。由于森林的年际变化不大且只有2015年的Landsat数据经过了校正，因此选用2015年同期的Landsat数据进行反照率反演。

3、反照率反演方法     (3.1) BRDF原型反演算法及WIDAS反照率的计算

AFX(Anisotropic Flat Index)是基于核驱动模型构造的角度指数，定义为白天空反照率和地表各向同性参数 f_iso的比，具体公式为(焦子锑 等，2011；Jiao 等，2012, 2014)

$ AFX = 1 + \frac{{{f_{\rm{vol}}}\left( \lambda \right)}}{{{f_{\rm{iso}}}\left( \lambda \right)}}{H_{\ker \_{\rm{vol}}}} + \frac{{{f_{\rm{geo}}}\left( \lambda \right)}}{{{f_{\rm{iso}}}\left( \lambda \right)}}{H_{\ker \_{\rm{geo}}}} $

(2)

式中，H_{ker_vol}，H_{ker_geo}为体散射核和几何光学核的双半球积分值。从式(2)可以看出，AFX是被各向同性参数标准化的几何光学权重参数和体散射权重参数与两个核的积分的线性表达式。

研究表明AFX较好地表征地表基本散射类型的变化，可由核驱动模型3个参数 f_iso，f_vol，f_geo直接求得AFX，并以此作为分类指标建立BRDF原型，用来表征地表BRDF形状的变化，具体方法可参见Jiao等人(2014)的文献。

本文参考Jiao等人(2014)提出的BRDF原型的建立方法，针对根河区域的地表各向异性变化特点，采用MODIS的BRDF参数(MCD43A1)数据，首先提取了相应的BRDF原型，具体参数信息如表2所示，其中F_iso′、F_vol′和F_geo′为标准化的模型参数。

基于核驱动模型，BRDF原型作为先验知识的反演算法。反演过程如下：根据核驱动模型前向模拟得到与观测数据(ρ_i)观测几何相同的方向反射率(ρ_i′)，采用最小二乘拟合法，最小化拟合误差，从而得到调整系数a

$ a = \frac{{\displaystyle\sum\limits_{j = 1}^n {{\rho _j} \cdot \rho _j'} }}{{\displaystyle\sum\limits_{j = 1}^n {{{\left( {\rho _j'} \right)}^2}} }} $

(3)

这样, 观测数据的BRDF就可以表示为

$ {\rm{BRDF}} = a \cdot {\rm{BRDF}}' $

(4)

拟合误差RMSE可以表示为

$ {\rm{RMSE}} = \sqrt {\frac{1}{{n - 1}}\sum\limits_{k = 1}^n {{{\left( {{\rho _k} - a\rho _k'} \right)}^2}} } $

(5)

得到RMSE最小的BRDF原型作为数据的先验BRDF；然后，通过式(3)确定调整系数a；最后，黑天空BSA和白天空WSA反照率与选定的最优原型的黑天空BSA′和白天空WSA′反照率之间的关系可以表示为

$ {\rm{BSA}}\left( {{\theta _i},{\textit{λ}} } \right) = a \cdot {\rm{BSA}}'\left( {{\theta _i},{\textit{λ}}} \right) = a \cdot \sum\limits_k {{f_k}{h_k}\left( {{\theta _i}} \right)} $

(6)

$ {\rm{WSA}}\left( {\textit{λ}} \right) = a \cdot {\rm{WSA}}'\left( {\textit{λ}} \right) = a \cdot \sum\limits_k {\left( {{f_k}\left( {\textit{λ}} \right){H_k}} \right)} $

(7)

式中，f_k是3个核参数，h_k(θ_i)是相应的核在观测半球的积分值，它们是太阳天顶角的函数；H_k是核在入射和观测半球的积分值。RossThickChen和LiSparseR对应的在双半球空间的积分值分别是0.189184和−1. 377622。

将每个波段的4种BRDF原型分别拟合WIDAS相应波段的多角度数据，将拟合误差RMSE最小的BRDF原型作为先验知识，即可得到WIDAS多角度数据对应波段的光谱反照率。最后，得到的波段反射率基于Liang(2001)提出的波段转换系数进行合成，即可得到相应的宽波段反照率数据。由于根河区域WIDAS多角度数据目前只有绿、红和近红波段数据可用，所以宽波段反照率转换系数采用AVHRR的二阶多项式转换公式(Liang，2001)，其表达式为

$ \begin{aligned} {\alpha _{{\rm{short}}}} = & - 0.3376\alpha _1^2 - 0.2707\alpha _2^2 + 0.7074{\alpha _1}{\alpha _2} + \\ & 0.2915{\alpha _1} + 0.5256{\alpha _2} + 0.0035 \end{aligned} $

(8)

式中，α₁、α₂为红和近红外波段的光谱反照率。

    (3.2) Landsat反照率的计算

Landsat数据具有长时间序列的全球观测，因此很多研究者基于这一数据发展高分反照率算法。本研究使用Landsat 8经过大气校正和地形辐射校正的地表反射率产品数据(Hu 等，2014)。

本文采用Shuai等人(2011)的算法，对Landsat数据的反照率进行计算，由于根河试验区域土地覆盖类型以针叶林为主，具有较好的地表一致性，因此本文不再进行Landsat地表反射率数据分类。首先将对应根河飞行区域30 m分辨率同时相8月15日的Landsat数据(63×91，2831个像元)的UTM投影转换成MODIS数据的Sinusoidal(SIN)投影，再聚合到500 m。然后利用与Landsat数据时空对应的高质量(QA=0)MODIS BRDF数据MCD43A1(4×6，24个像元)计算AN(Albedo-to-Nadir)系数，过程中假设MODIS在Landsat分辨率下是纯像元，则Landsat的AN系数就近似为MODIS数据的AN系数。最后，AN系数乘以Landsat反射率就得到30 m分辨率Landsat光谱反照率，再转换为与WIDAS数据相同的UTM投影，便于WIDAS反照率对其进行验证。Landsat窄波段到宽波段的转换系数得到地表短波宽波段反照率的公式(Liang，2001)

$ \begin{aligned} {\alpha _{{\rm{short}}}} =& 0.356\alpha _1^2 + 0.130{\alpha _3} + 0.373{\alpha _4} + \\ & 0.085{\alpha _5} + 0.072{\alpha _7} - 0.0018 \end{aligned} $

(9)

计算过程如图3所示。

4、分析与验证

不同于黑河区域的农作物试验区，根河森林试验区内由于地形和树高等影响，通量塔站点和临时反照率表的设置较为困难。因此，其数量相对有限。同时，根河区域WIDAS没有多时相的飞行数据，无法进行站点的时间序列的验证。在2008年和2012年黑河农田区域，已有研究发展了针对WIDAS传感器的反照率反演算法，并进一步验证算法的鲁棒性(Jiao 等，2012, 2014；张虎 等，2013)。

试验区域内的一个地面通量塔验证点，其周围的地表类型为针叶林，地形较为平坦。虽然一个能量塔的验证数据具有随机性，但仍作为一定参考。区域内通量塔观测的反照率为0.177，Landsat计算的反照率为0.150，原型反演结果为0.168。很明显，原型反演方法的结果与真实观测更为接近，但Landsat反演的反照率相对于通量塔的测量结果有些低估，平均低估0.027，这个结果与估算的两种遥感数据的均方根误差(即RMSE为0.02)较为接近。因此，该通量塔的结果表明，BRDF原型算法反照率反演结果相对较为准确，并且已在黑河飞行区域经不同时期同步地表反照率表测量数据的验证，不确定性更小，在该研究中可以被近似认为真值。

将BRDF原型反演算法计算的1 m空间分辨率的WIDAS反照率(图4)，重采样到30 m(图5)。上述3.2节方法得到的根河飞行区域30 m分辨率的Landsat反照率结果(图6)。图中绿色五角星是通量塔的位置。很显然，1 m WIDAS地表反照率(图4)提供了非常丰富的纹理与细节。目视比较图5和图6，可以看出，BRDF原型算法得到的地表反照率(图5)和Landsat反照率结果(图6)整体的趋势较为一致，相对于WIDAS反照率的反演结果，Landsat反照率总体略有偏低。由图4可以清晰地看到，图中右下角有条带状较低的反照率值，Landsat反照率同样存在这种情况。这是因为该区域存在与条带走向垂直的坡度变化，且飞行时的太阳天顶角为40°，方位角为50°左右。Landsat数据虽对地形起伏明显的区域(坡度大于8°)进行了校正(Wang 等，2016)，但仍有部分区域地形的起伏产生阴影导致反照率较低。

二者的差异可能存在多方面的原因。Landsat反照率在计算过程中更可能引入随机误差和系统误差，随机误差包括在特定气溶胶条件下，Landsat反射率可能存在较小的过校正情况(Masek 等，2006)。系统误差可能在于窄波段向宽波段转换，以及500 m的MODIS BRDF参数应用到30 m的Landsat数据过程中存在尺度效应问题。另外，MODIS和Landsat传感器的带通差异也会引进一定的误差(Shuai 等，2011)。同时，对于单个角度的Landsat观测，对于可能的误差影响(如，亚像元云的影响)，该算法会将这种误差完全传递到反照率结果中(Wanner 等，1995)。多角度WIDAS机载数据在一定程度上会较好避免这种误差的传递(Jiao 等，2015)。

由于飞行区域比较有限，为了更深入地分析根河试验区整体的反照率情况，选取了较大区域50°53′—50°58′N，121°26′—121°35′E，包含318×347个Landsat像元，基于Shuai等人(2011)的算法计算反照率，如图7所示，图中绿色多边形即为上述验证的飞行区域，绿色圆点代表通量塔的位置。整个区域内地表类型是针叶林，有一定的地形起伏，飞行区域相对整个区域较为均匀且地形起伏较小，反照率值连续性较好。区域内的反照率值有90%集中在[0.1, 0.35]，这与飞行区域的反照率范围基本一致，其中小部分较大或较小值主要是由于地形和云等的影响。

为了更直观地对反演结果进行比较，图8进一步定量比较了Landsat反照率和BRDF原型算法反演的反照率散点图。图8中颜色分别表示在边长为0.005的方格内，将落在方格内的反照率的个数从高到低五等分并赋为红、黄、绿、蓝和灰色的密度散点。

BRDF原型算法反演的反照率和Landsat反照率相比(图8)，偏差基本在0.05以内，散点图集中。表明应用Shuai的算法反演的Landsat反照率结果与BRDF原型算法的结果在飞行区域总体上一致性较好。均方根误差为0.02，Bias为0.0057，但Landsat反照率结果相对于WIDAS机载的反照率结果，总体上有些低估。主要原因如上述分析所示，可能是两种遥感数据的一些可能的系统误差，如大气校正，传感器的辐射定标与校正的差异，会很大程度上影响到比较结果。

5、结　论

地表反照率的验证是遥感参数验证的核心内容。目前地表反照率主要验证方法的发展主要是基于通量塔观测的地表反照率数据，通过分析这些数据的代表性，然后与星载数据反演的地表反照率进行直接比较(Román 等，2009)。然而，机载飞行数据可以在短时间内获取面源的多角度数据，为星载反照率的验证提供重要手段。机载WIDAS多角度飞行试验，在黑河和根河试验中，历经近十年的改进与发展，针对不同的地面目标(从农作物到森林)，发展和验证了针对小观测角度的地表反照率反演算法。根河试验中进一步探讨了BRDF原型反演算法应用于森林试验区验证Landsat反照率的情况。

本研究以核驱动模型为基础，应用BRDF原型算法对根河地区的反照率进行反演。然后，进一步将反演结果对该区域Landsat反照率进行验证。结果表明，在观测角度最大仅为25°的情况下，BRDF原型算法得到的WIDAS机载反照率结果与Landsat反照率进行比较，总体均方根误差为0.02，偏差为0.0057，满足对反照率绝对精度的要求。

本文的工作主要体现在以下几个方面：(1)针对根河实验区域森林地表覆盖类型，提取了与之对应的MODIS BRDF原型作为先验知识，因此，更能代表该区域林地的BRDF特征；(2)应用机载WIDAS数据对Landsat反照率进行初步验证，丰富了Landsat卫星反照率数据的验证形式。

由于目前根河区域的数据只有部分可用，在以后的工作中考虑收集更多的数据和验证点数据。另外，如何减小Landsat反照率在计算过程中的随机误差和系统误差，解决MODIS BRDF参数产品和Landsat以及WIDAS数据之间的尺度效应的影响，也是一个有待解决的科学问题。

志 谢 感谢WIDAS飞行团队的刘强、肖青、闻建光、曾也鲁、曹彪等各位老师提供数据和进行数据的预处理工作。Landsat校正数据是刘良云老师提供。在此一并感谢！