基于HJ-1A/B CCD地表反照率估算方法比较与验证

樊宪磊， 阎宏波， 瞿瑛

(东北师范大学地理科学学院，长春 130024)

0 引言

地表反照率是辐射与能量平衡的重要特征参量，被广泛应用于中长期天气预报和全球气候变化研究[1]。目前，已经基于卫星遥感数据生成了大量反照率数据集[2]，如MODIS[3]，GLASS[4-5]和CLARA[6]等。但是，这些反照率数据集大多基于中低空间分辨率卫星观测数据生成，因此空间分辨率相对较低，难以准确地描述景观破碎度较大区域的地表反照率时空变化情况[2]。而高空间分辨率卫星遥感数据可以提供丰富的地表特征信息，使得基于高空间分辨率卫星遥感数据的地表反照率估算方法[7-8]成为研究中新的增长点[5]。

环境减灾小卫星(HJ-1A/B)可以提供大幅宽、短重访周期的30 m空间分辨率CCD影像数据，为生成高空间分辨率地表反照率数据集提供了可靠数据源。孙长奎等[9]基于POLDER BRDF数据集，通过分地类和分格网回归的方式，建立了HJ-1A/B CCD地表反射率与宽波段反照率之间的回归关系； He等[10]提出了一种基于HJ-1卫星大气层顶反射率(top of atmosphere reflectance,TOA)数据的反照率直接反演算法(direct estimation algorithm，DEA)，证明了基于HJ-1A/B CCD TOA进行反照率直接估算的可行性； Gao等[11]在Shuai等[7]Landsat反照率产品算法的基础上，提出了一种将HJ-1A/B卫星地表反射率(surface reflectance,SUR)数据与MODIS双向反射分布函数(bidirectional reflectance distribution function,BRDF)产品相结合的反照率估算方法，通过建立HJ-1A/B CCD和MODIS SUR比值关系来生成高空间分辨率的地表反照率数据集； 张虎等[12]将MODIS核系数(MODIS kernel coefficients,MKC)产品作为先验知识，基于贝叶斯推论对HJ-1 CCD地表反照率进行了反演，研究结果表明该方法能够改善其他方法中出现的马赛克现象。

这些方法都能够实现基于HJ-1A/B CCD数据生成高空间分辨率的反照率数据集，但其图像精细度和估算精度还缺乏系统性的评价和比较验证。本研究从图像精细度和估算精度2个方面，对2种基于HJ-1A/B CCD数据的地表反照率遥感估算方法,即DEA和基于MKC的估算方法，进行了比较和验证分析，为地表反照率数据集生成方法选择和算法优化改进研究提供参考。

1 地表反照率估算方法

1.1 DEA

DEA是基于BRDF数据库建立TOA/SUR与宽波段反照率之间的回归关系[5]，从而实现基于单一角度观测数据直接估算地表宽波段反照率的方法[13]。该方法首先基于POLDER BRDF数据集通过插值和筛选构建SUR训练数据集； 然后，再通过波段转换、角度积分、窄波段向宽波段转换和大气辐射传输模型等过程，构建包含TOA/SUR和对应地表宽波段反照率的训练数据集； 最后，通过分地类(分为植被、土壤和冰雪3类)和分角度格网(太阳天顶角和观测天顶角2°间隔，相对方位角5°间隔)回归的方法构建TOA/SUR与宽波段反照率之间的线性回归关系，生成DEA查找表。根据太阳/观测几何信息以及SUR数据查找对应的直接反演系数，可以直接估算得到地表宽波段反照率。该方法不受MODIS反照率产品生成算法中地表BRDF特性需要在较短时间窗口内(如16 d)不发生显著变化假定条件的限制，因此更适用于监测地表BRDF特性发生快速变化时(如降雪、融雪、森林火灾和农作物收割等)地表反照率的时空动态变化过程[5]。DEA具体流程如图1所示。

图1 DEA流程Fig.1 Flowchart of the DEA

DEA有2种形式，可以通过SUR或者TOA来估算宽波段反照率，分别记做DEA-SUR和DEA-TOA。其中DEA-SUR算法可以基于SUR直接估算地表宽波段反照率，即

(1)

(2)

式中：αbs和αws分别为地表宽波段的黑空和白空反照率；ρi(θs,θv,φ)为第i波段的SUR；θs,θv和φ分别为太阳天顶角、观测天顶角和相对方位角；θbs为局地正午太阳天顶角；ci和c0分别为DEA查找表中的回归系数。

DEA-TOA算法可以基于TOA直接估算地表宽波段反照率，即

(3)

(4)

通过比较可以发现DEA-SUR和DEA-TOA具有相似的估算精度，为了降低比较分析的复杂性，仅将DEA-SUR作为本研究比较分析评价的算法。

1.2 基于MKC的地表反照率估算方法

该方法假定MODIS纯像元与所包含的高空间分辨率亚像元具有近似的BRDF性质，因此可以基于MKC估算高空间分辨率的地表反照率[7]，估算方法流程如图2所示。

图2 基于MKC的估算方法流程Fig.2 Flowchart of the method based on MKC

首先通过对影像数据进行图像分类，区分出纯像元和混合像元。对于纯像元，假定低空间分辨率的MODIS像元与内部的高空间分辨率HJ-1A/B CCD像元具有相同的BRDF性质，以MODIS BRDF产品(MKC和MCD43A1产品)作为先验知识，通过MODIS的反照率与反射率比值等于内部的HJ-1A/B CCD亚像元反照率与反射率的比值关系估算得到高空间分辨率的地表反照率，即

(5)

(6)

式中：αbs,M(θbs)为局地正午太阳天顶角为θbs时的MODIS黑空反照率；αws,M为MODIS白空反照率；λ为波长；Rλ,M(Ω)为入射/观测几何为Ω时的MODIS波段反射率；αbs,HJ(θbs)为局地正午太阳天顶角为θbs时的HJ-1A/B CCD黑空反照率；αws,HJ为HJ-1A/B CCD的白空反照率；Rλ,HJ(Ω)为观测几何为Ω时的HJ-1A/B CCD波段反射率；Cλ为反照率与波长为λ的波段反射率的比值。

对于混合像元，MODIS的BRDF代表该像元内各种地表覆盖类型的BRDF加权组合，参照周围具有相同地表覆盖类型的MODIS纯像元BRDF，构造3 km×3 km窗口内的基于HJ-1A/B地表类型数据与地表反照率的经验关系，通过拟合获得加权经验系数，最终得到高空间分辨率的地表反照率，即

αλ,HJ=CoRλ,HJ(Ω)，

(7)

式中αλ,HJ为HJ-1A/B CCD的地表反照率；Co为拟合得到的经验系数。最终宽波段的地表反照率可以通过窄波段向宽波段转换公式计算得到[11]。

2 比较与验证评价方法

2.1 图像精细度评价方法

本研究采用目视判读和图像清晰度指数评价基于HJ-1A/B CCD数据的地表反照率遥感估算方法的图像精细度，并采用MODIS数据作为图像精细度评价的参考数据。通过目视判读方式对比DEA-SUR，MKC算法和MODIS地表反照率产品的图像精细度和地物识别能力。为了更加客观地评价基于HJ-1A/B CCD数据估算得到的高空间分辨率地表反照率数据对空间变化的刻画能力，还计算了DEA-SUR，MKC算法与MODIS反照率产品的图像清晰度指数。采用拉普拉斯梯度函数[14]进行图像边缘检测，将图像中检测为边缘的拉普拉斯算子卷积的均值作为衡量图像清晰度的指标，值越大表明图像细节越丰富，清晰度越高。在拉普拉斯梯度函数中定义的拉普拉斯算子L为

(8)

图像清晰度D(f)定义为

(9)

式中：G(x,y)为像素点(x,y)处拉普拉斯算子的卷积；T为边缘检测的阈值(本研究设T=0.01)；n为待评价图像像元总数。

2.2 基于站点观测数据的验证方法

为了评价和比较基于HJ-1A/B CCD数据估算地表反照率算法的精度，使用US-MMS、长岭(CN-Cng)、盈科(Yingke)和纳木错(Namco) 4个站点(表1)观测数据对DEA-SUR和MKC算法进行了比较验证，并使用MODIS地表反照率产品作为参考数据。

表1 地表反照率地面观测站点信息Tab.1 Information of the in situ measurements sites

其中US-MMS和长岭站点为Fluxnet全球通量观测站点，盈科和纳木错站点数据为中国科学院寒区旱区环境与工程研究所和青藏高原研究所实验观测数据。这些站点的地表覆盖类型包括落叶阔叶林、草地、农田和高寒草甸，观测站点的Google Earth影像如图3所示。考虑到地面站点观测与卫星遥感影像观测的空间尺度存在较大的差异，对4个站点地表反照率随空间分辨率的变化也进行了评价分析。

(a) US-MMS (b) 长岭 (c) 盈科 (d) 纳木错

图3 地表反照率观测站点的Google Earth影像示意图

Fig.3GoogleEarthimageriesinsitumeasurementssites

选取站点每天局地正午时刻观测的太阳上行辐射和下行辐射的比值作为地表反照率观测值(蓝空反照率)。基于卫星观测数据估算得到的黑空和白空反照率需要基于天空散射比进行线性加权计算得到蓝空反照率[15]，即

αblue(θbs)≈[1-S(τ,λ)]αbs(θbs)+S(τ,λ)αws，

(10)

式中：αblue(θbs)为蓝空反照率；S(τ,λ)为大气散射光比例；τ为气溶胶光学厚度。在本研究中，大气散射光比例通过基于6S程序建立的查找表计算获得。为了客观评价和验证不同算法的精度，绘制了站点观测值与基于HJ-1A/B CCD数据估算的蓝空反照率的时间序列图和散点图，计算均方根误差(root mean squared error，RMSE)、偏差值(Bias)和决定系数(R2)来评价和验证DEA-SUR和MKC算法的估算精度。

3 结果与分析

3.1 图像精细度评价结果

环境减灾小卫星星座包括HJ-1A，B和C星3颗卫星，单颗卫星的重访周期为4 d，当HJ-1A和HJ-1B这2颗卫星组网后重访周期可达到2～4 d。其中HJ-1A/B每颗卫星上分别搭载了2台多光谱CCD传感器，分别都包含蓝光、绿光、红光和近红外4个波段。为了评价基于HJ-1A/B CCD影像生成的地表反照率的图像精细度，选择了2009年8月14日在张掖市获取的无云覆盖、成像质量良好的HJ-1A CCD影像(图4)作为图像精细度评价数据，下载自中国资源卫星应用中心陆地观测卫星数据服务平台。

图4 张掖市HJ-1A CCD B4(R)，B3(G),B2(B)假彩色合成影像Fig.4 False color imagery with B4(R),B3(G),B2(B)of HJ-1A CCD in Zhangye City

本研究使用MCD43A3产品作为图像精细度和估算精度评价的参考，该产品的时间分辨率为1 d，空间分辨率为500 m，产品共包含MODIS 1—7波段和3个宽波段(0.3～0.7 μm(可见光)、0.7～5.0 μm(近红外)、0.3～5.0 μm(短波红外))的白空和黑空反照率(局地正午)。MODIS地表反照率产品是美国国家航空航天局 MODIS团队使用陆地表面各向异性二向性反射模型算法(algorithm for MODIS bidirectional reflectance anisotropy of the land surface，AMBRALS)[16]生成的。AMBRALS算法假定地表BRDF在一个时间窗口(如16 d)内保持不变，通过累积16 d内无云覆盖的多角度、多波段的MODIS观测数据，通过大气校正、BRDF角度建模和窄波段向宽波段转换等步骤生成MODIS地表反照率产品[3]。

基于HJ-1A/B CCD数据的地表反照率估算方法(DEA-SUR和MKC)与MODIS地表反照率产品在张掖市实验区的比较如图5所示。

(a) DEA-SUR地表白空反照率 (b) MKC地表白空反照率 (c) MODIS地表反照率产品

(d) 图(a)矩形区域放大 (e) 图(b)矩形区域放大 (f) 图(c)矩形区域放大

图5 张掖市实验区地表反照率数据集(白空反照率)比较

Fig.5Comparisonofthesurfacealbedodatasets(white-skyalbedo)inZhangyeCity

通过目视判读对比可以发现，DEA-SUR和MKC算法估算结果在图像精细度方面要显著优于MODIS地表反照率产品，比较结果表明基于HJ-1A/B CCD数据的地表反照率估算方法能够清晰地刻画出空间破碎度较大区域的地表反照率空间变化情况。基于清晰度指数评价结果，MKC算法反演的结果检测到的边缘最丰富(清晰度指数为0.192)，但是内部存在噪声点，由于使用了500 m空间分辨率的MKC产品，在一些区域存在明显的马赛克现象(图5(e))； DEA-SUR算法估算结果的细节丰富度次之(清晰度指数为0.120)，但是检测到的边缘清晰，细节丰富，噪声点相对较少，无MKC算法中存在的马赛克现象； MODIS产品的空间分辨率相对较低，只能检测出图像的大致轮廓信息(清晰度指数为0.094)，图像精细度最差。

3.2 基于站点观测数据的验证结果

地面观测站点的地表反照率随空间分辨率(黑色矩形框代表了10种不同的空间分辨率)的变化如图6所示，图中矩形从内向外依次代表30，90，150，210，270，330，390，450，510和1 050 m的10种空间分辨率像元覆盖区域。图7统计了站点在不同空间分辨率的地表反照率均值，其中US-MMS和盈科站点统计时段为植被生长季(5—10月)，纳木错站点的统计时段为冬季(12月—次年2月)，长岭站点统计时段为全年。使用不同空间分辨率的地表反照率估算值与站点实测值计算RMSE作为衡量站点空间代表性的指标，图7中误差线长度即代表RMSE大小，表征了不同空间分辨率像元反照率估算值与实测值的变异程度，反映了站点观测值所能表征的空间尺度。统计结果表明US-MMS、长岭、盈科和纳木错站点在无雪时段都具有较好的空间代表性(US-MMS站点RMSE为0.03左右，长岭站点RMSE为0.03左右，盈科站点RMSE为0.02左右，纳木错站点RMSE为0.04左右)，但长岭站点在积雪覆盖时段空间代表性稍差(长岭站点积雪覆盖时段RMSE为0.08左右)。

(a) US-MMS (b) 长岭 (c) 盈科 (d) 纳木错

图6 地表反照率随空间分辨率的变化

Fig.6Variationofsurfacealbedowiththespatialresolution

图7 地表反照率随空间分辨率的变化统计Fig.7 Variations of surface albedo with spatial resolutions

在地表覆盖类型为落叶阔叶林的US-MMS站点，地表反照率在0.1～0.2之间变化。由于HJ-1A/B在国外有效观测数据数量较少，因此从2008—2014年间，除去有云覆盖的影像外符合要求的数据只有15 d。基于有限的观测数据验证结果表明(图8)，DEA-SUR与地面站点观测数据的一致性较好(RMSE为0.015)，MKC算法与地面站点观测数据的一致性次之(RMSE为0.018)。由于地表下垫面类型为落叶阔叶林，反照率表架设在较高的通量观测塔上，因此该站点代表性和空间均一性较好，DEA-SUR，MKC和MODIS反照率产品均与站点观测数据具有较好的一致性，可以证明DEA-SUR和MKC算法都能够有效地估算地表反照率。

(a) 时间序列 (b) 散点图

图8 DEA-SUR，MKC和MODIS反照率与US-MMS站点观测值比较

Fig.8ComparisonofthealbedobyDEA-SUR,MKCandMODISproductwiththeinsitumeasurementsatUS-MMSsite

在地表覆盖类型为草地的长岭站点，地表反照率在夏季的最低值在0.15左右，秋季草枯黄后反照率逐渐升高，冬季降雪时地表反照率可升高至0.5～0.6(图9)。在无积雪覆盖的时段，DEA-SUR和MKC算法与地面站点观测数据具有较高的一致性(RMSE均为0.023)； 而在积雪覆盖时段，2种算法和MODIS产品均有不同程度的高估(DEA-SUR的Bias为0.141，MKC的Bias为0.130，MODIS产品的Bias为0.103)。通过Google Earth影像目视判读，这主要由于站点积雪被人为踩踏影响，使得地表反照率的空间变化更为剧烈，站点的空间代表性变差，因此站点观测结果与卫星观测数据不具备可比性。

(a) 时间序列 (b) 散点图

图9 DEA-SUR，MKC和MODIS反照率与长岭站点观测值比较

Fig.9ComparisonofthealbedobyDEA-SUR,MKCandMODISproductwiththeinsitumeasurementsatCN-Cngsite

在地表覆盖类型为农田的盈科站点，地表反照率在2009年5—8月间反照率变化幅度不大(0.15～0.2)(图10)。MODIS产品与地面站点观测结果一致性最好(RMSE为0.007)，反照率变化趋势较为稳定和平滑。而MKC算法(RMSE为0.023)和DEA-SUR算法(RMSE为0.027)的精度相当，估算误差相对MODIS地表反照率产品较大。这主要是由于MODIS产品是基于16 d观测数据反演得到的，相当于对反照率数据集进行了时间序列滤波和平滑处理。

(a) 时间序列 (b) 散点图

图10 DEA-SUR，MKC和MODIS反照率与盈科站点观测值比较

Fig.10ComparisonofthealbedobyDEA-SUR,MKCandMODISproductwiththeinsitumeasurementsatYingkesite

在地表覆盖类型为高寒草甸的纳木错站点，在无积雪覆盖时地表反照率约为0.2，在冬季降雪时地表反照率可升高至0.8左右，在积雪融化时反照率快速降低(图11)。在2009年12月—2010年2月，2种基于HJ-1A/B CCD数据的地表反照率估算方法精度相当，在积雪覆盖时段，DEA-SUR和MKC算法的RMSE均为0.058，而MODIS地表反照率产品的估算误差相对较大(RMSE为0.134)。在无积雪覆盖时，DEA-SUR和MKC算法的RMSE分别为0.032和0.041。

(a) 时间序列 (b) 散点图

图11 DEA-SUR，MKC和MODIS反照率与纳木错站点观测值比较

Fig.11ComparisonofthealbedobyDEA-SUR,MKCandMODISproductwiththeinsitumeasurementsatNamcosite

4 结论与讨论

本研究从图像精细度和估算精度2个方面对基于HJ-1A/B CCD数据的2种反照率估算方法(DEA和MKC)进行了比较与验证分析。通过目视判读和清晰度指数评价了高空间分辨率反照率产品的图像精细度，基于US-MMS、长岭、盈科和纳木错4个地面观测站点数据验证了2种地表反照率估算方法的精度。研究得出的主要结论如下：

1)基于目视判读和清晰度指数的图像精细度评价分析结果表明，基于HJ-1A/B CCD数据的反照率估算结果相比于MODIS反照率产品能够反映出更为丰富的反照率空间分布与变化细节，其中MKC和DEA-SUR的图像清晰度相当，但是MKC算法生成的反照率数据集存在明显的马赛克现象，而DEA-SUR算法则能够较好地解决该问题。

2)基于站点观测数据的比较和验证分析结果表明，在无积雪覆盖时，采用DEA-SUR算法在4个站点的RMSE： US-MMS为0.015，长岭为0.023，盈科为0.027，纳木错为0.032； 采用MKC算法在4个站点的RMSE： US-MMS为0.018，长岭为0.023，盈科为0.023，纳木错为0.041，2种方法与地面站点观测数据的一致性相当。而在积雪覆盖时，2种方法的估算精度与无积雪覆盖时相比都存在明显下降，DEA-SUR在长岭和纳木错站点的RMSE分别为0.155和0.058，MKC算法在长岭和纳木错站点的RMSE分别为0.135和0.058。

3)基于HJ-1A/B CCD数据估算的地表反照率数据存在更多的抖动和噪声，估算精度总体略低于MODIS产品。这主要是由3方面原因造成的： ①MODIS反照率产品是基于16 d卫星观测数据生成的，相当于对反照率进行了时间滤波和平滑处理，因此相比DEA-SUR和MKC算法更为稳定和平滑； ②由于HJ-1A/B CCD只有可见光—近红外4个波段，缺少短波红外波段，因此总体估算精度会略低于具有可见光—短波红外7个波段的MODIS数据的估算精度； ③因为DEA-SUR和MKC算法是基于低空间分辨率BRDF先验知识(POLDER BRDF数据集和MCD43A1产品)，而研究表明BRDF信息随空间分辨率存在着显著的变化，因此低空间分辨率的BRDF先验知识也会引入估算误差。因此在未来需要针对上述3方面问题寻找改进方案，不断提高基于HJ-1A/B CCD数据生成的地表反照率的图像精细度和估算精度。