陈辉,杨啸,李营,张丽娟,王中挺
(1.生态环境部卫星环境应用中心,北京 100094;2.中国人民大学 环境学院,北京 100872)
高分一号(GF-1)卫星作为中国高分系列的首发星,自2013年4月成功发射以来,成为对地观测的主要卫星数据源之一,在我国国土、农业、环境等遥感监测中得到了广泛的应用。随着卫星遥感技术的不断发展和遥感监测指标定量化需求的日益增强,定量遥感逐渐成为遥感应用发展的主要方向,然而由于大气效应的存在影响了遥感图像的可用性,需要进行大气校正获得地表真实反射率以提高产品质量[1]。因此,准确的大气校正是卫星遥感定量化应用的前提与基础。
目前,国内外已提出多种大气校正方法模型,主要包括图像特征模型、地面线性回归经验模型、大气辐射传输理论模型等,其中大气辐射传输理论模型以数学和物理理论为基础,较为科学合理地呈现了太阳辐射在太阳-地球-卫星系统中的传输过程,物理意义较为明确且精度相对较高,为目前使用最为广泛的大气校正方法模型[2]。Peng等[3]采用6S(second simulation of the satellite signal in the solar spectrum)模型、郑盛等[4]采用MODTRAN模型、刘伟刚等[5]采用FLAASH模型对TM、CBERS、HJ等卫星进行大气校正研究,结果表明,利用辐射传输模型进行大气校正可有效减少或消除大气气溶胶对卫星影像识别地物的干扰。近年来,刘佳等[6]针对GF-1卫星的数据特点采用6S模型进行大气校正算法研究,并与FLASSH校正结果进行比较,结果表明,6S模型在大气校正方面具有高精度、高效率的计算优势。
大气辐射传输模型法已经成为当前主流的遥感影像大气校正方法之一。这种方法有一个重要前提,即大气参数的获取,其中大气气溶胶光学厚度(aerosol optical depth,AOD)是最主要的大气参数之一[7]。上述大部分大气校正研究中,一般采用同一个能见度或者AOD对整幅影像进行辐射传输计算实现大气校正,这从一定程度上忽略了区域上大气环境的空间变化,容易造成“清洁”地区反射率大气校正过多而“污染”地区反射率大气校正不足的现象。GF-1卫星搭载了4个宽视场传感器(wide field of view,WFV),每个相机的视场角约16°,幅宽约为200 km,在大气校正时,不仅要考虑到大气环境情况的时间变化,还应考虑到其空间上的差异。文献[8-9]结合MOD09地表反射率,利用GF-1的WFV和PMS相机进行AOD反演并计算地表反射率,但是由于缺少红外波段的支持,导致AOD结果精度有限,从而影响了地表反射率的计算精度。
本文基于GF-1卫星WFV相机数据特征,从MOD04产品中提取大气校正的关键输入参数,并采用6SV矢量辐射传输模型模拟计算大气参数,建立基于GF-1卫星WFV数据的逐像元大气校正方法及技术流程。
作为我国高分辨率对地观测系统系列的首发卫星,GF-1卫星搭载了4个WFV相机,具有高空间分辨率(星下点16 m×24 m)和较大的幅宽(4 km×200 km),重访周期为4 d[10]。每台WFV相机都能够获取4波段多光谱卫星遥感影像,其波段响应曲线和设计参数如图1和表1所示[11]。
图1 GF-1号WFV1相机4个波段的响应函数
表1 GF-1卫星 WFV相机设计参数
6S模型是美国马里兰大学Vermote等[12]在卫星信号在太阳光谱中的模拟(simulation of the satellite signal in the solar spectrum,5S)模型基础上,发展用于模拟0.25~4.0 μm间无云大气条件下电磁辐射在太阳-目标物-传感器系统中传输过程变化的模型。在6S模型的基础上,研究人员考虑了大气和地表的极化特性,进一步发展了矢量辐射传输模型的卫星信号在太阳光谱中的二次模拟矢量模型(second simulation of the satellite signal in the solar spectrum vector version,6SV),并对散射角度、波长节点、气溶胶垂直廓线等特征参数进一步丰富和完善。用户可根据实际需要对上述参数进行自定义设置,提高了辐射传输模型的灵活性、合理性和实用性。6SV通过设定输入特定的大气、光谱和几何观测条件,计算气体吸收、气溶胶和大气分子散射对地表-大气耦合系统的影响,从而模拟出地表对太阳辐射情况[13]。本研究采用2015年6月发布的6SV2.1版本软件包。
为对比分析6S和6SV对GF-1卫星WFV相机的大气校正效果,设定地表类型(植被)和太阳及卫星观测几何信息(太阳天顶角45°,观测天顶角为0°),利用6S和6SV模型分别模拟了不同AOD条件下WFV相机4个波段的表观反射率(图2)。从模拟结果可以看出,当AOD在1以内时,GF-1卫星WFV相机4个波段6S和6SV 2种辐射传输模型模拟的表观反射率较为接近,随着AOD的升高,2种辐射传输模型模拟的表观反射率开始出现差异。尤其当AOD大于2时,6S模拟的表观反射率开始出现“梯度变化”。例如,当AOD在1.8~2.3和2.6~3.2之间变化时,band1的表观反射率基本不变;当AOD在2.6~3.2和3.5~4之间变化时,band2的表观反射率变化幅度非常小;当AOD在2.9~3.2和3.5~4.5之间变化时,band3的表观反射率变化幅度非常小;当AOD在2~5之间变化时,band4的表观反射率变化幅度非常小;而6SV模拟的表观反射率则仍表现出连续性变化。这说明2种辐射传输模型在低AOD情况下对GF-1卫星WFV相机的模拟能力相当,但在AOD较高时,6S对AOD变化反应较为迟钝,6SV的模拟结果对AOD变化更为灵敏,这说明6SV进行大气校正时对不同AOD反应更精确。因此,选用6SV作为GF-1卫星WFV相机大气校正的辐射传输模型。
图2 基于GF-1卫星WFV相机参数的6S和6SV对AOD变化的敏感性模拟分析
美国国家航空航天局(National Aeronautics and Space Administration,NASA)基于搭载在Terra和Aqua卫星上的中分辨率光谱成像仪(MODIS)观测结果,提供全球尺度逐日的AOD数据。通过多年的算法改进和产品完善,目前最新版本的产品为2014年发布的C6版本,包括10 km和3 km 2种分辨率结果,即MOD04和MOD04_3k。C6版本较之前的C5版本在反演算法和产品质量上都有所提高,主要包括更高的覆盖率、更完善的瑞利光学厚度和气体吸收的假设、改进的植被指数、改进的质量控制方法以及更新的云掩膜等方面[14],并且大部分(约70%)的AOD数据结果误差在±(0.05+15%)以内[15]。NASA发布的MOD04产品为卫星遥感大气校正提供了丰富的气溶胶信息资料,其中MOD04_3k仅采用暗像元算法反演获取了3 km分辨率的AOD,MOD04则提供了包括暗像元、深蓝及2种算法联合生成的产品,其产品覆盖率和精度都达到较高的水平。因此,本研究选用MOD04产品,利用IDL设计程序从中提取暗像元-深蓝联合算法反演的AOD数据集(AOD_550_Dark_Target_Deep_Blue_Combined),并进行投影转换和拼接等预处理,为实现GF-1卫星WFV相机的大气校正提供参数。
假设地表为均匀朗伯体的情况下,卫星传感器接收到的表观反射率可用式(1)表示为大气和地表之间多次散射的结果。
(1)
式中:ρTOA为卫星接收到的表观反射率;θs与θv分别为太阳天顶角与观测天顶角;φ为太阳入射与卫星观测方向的相对方位角;ρ0为大气程辐射;S为大气整层的半球后向反射率;T(θs)为太阳到地面的散射透过率;T(θv)为地表到卫星传感器的散射透过率,一般T(θs)和T(θv)总是以乘积形式出现,因此,将T(θs)T(θv)可视为一个参数表示大气整层的散射透过率T;ρs为地表反射率。
而大气校正的主要目的是为了获取准确的地表反射率,根据式(1)进行分离变量可以获取地表反射率,计算如式(2)所示。
(2)
式中:地表反射率可表示为卫星表观反射率ρTOA、大气程辐射ρ0、大气散射透过率T和大气半球后向散射率S的函数。因此,要精确计算地表反射率就需要获取式(2)中的4个参数,其中,卫星表观反射率ρTOA可以通过辐射定标后的传感器信号计算获取,大气程辐射ρ0、大气散射透过率T和大气半球后向散射率S这3个参数实际上是大气状况的基本参数,主要与太阳-卫星观测几何及大气气溶胶相关。
本研究的主要思路即以区域AOD和GF-1卫星的观测几何参数为6SV模型输入参数,经过大气辐射传输模拟获取大气校正所需的参数,然后根据式(2)计算地表反射率。数据处理流程主要由辐射定标和表观反射率计算、AOD提取和校正、6SV模型模拟大气参数、计算地表反射率4个部分构成(图3)。
图3 数据处理流程图
根据GF-1卫星的过境时间选取辐射定标系数Gain和L0,将遥感影像像元亮度DN值根据式(3)转换为表观辐亮度L。
L=DN*Gain
(3)
根据从xml文件中获取的太阳天顶角θs,将表观辐亮度根据式(4)转为各波段的表观反射率。
(4)
式中:Eλ为大气层顶太阳辐照度;ds为日地距离修正因子。式(3)中所用的定标系数和式(4)中的太阳辐照度来自于中国资源卫星应用中心官方网站公布的GF-1卫星参数资料(http://218.247.138.119/CN/Downloads/dbcs/index.shtml)。
利用辐射传输模式结合GF-1卫星WFV相机的4波段光谱响应函数,根据观测几何(太阳天顶角、观测天顶角、相对方位角)和AOD,得到卫星遥感图像上不同区域的不同AOD下的大气参数(大气程辐射ρ0、大气层向下的半球反射率S和整层大气的透过率T),然后逐波段将获得的大气参数和WFV相机探测数据代入式(2),分别获取4个波段的地表反射率ρs,从而实现大气校正。
选取2015年5月25日过境华北地区的GF-1卫星WFV1相机观测资料(数据名称:GF1_WFV1_E115.8_N38.0_20150525_L1A0000826550)进行算法实验。为对大气校正结果进行对比分析,采用2种方法分别对同一景影像进行大气校正:①AOD均值大气校正,即计算卫星数据覆盖区域范围内的AOD平均值,利用平均AOD对整幅影像进行大气校正;②逐像元大气校正,即根据MOD04的区域AOD空间分布通过插值获取GF-1卫星WFV相机观测的每个像元的AOD值,并根据上述算法原理实现对GF-1影像逐像元大气校正。之后,将2种大气校正结果进行对比分析。
从实验结果可以看出,2015年5月25日,GF-1卫星WFV相机在华北地区观测到的影像(图4(a))整体偏暗且呈“雾”状模糊效果,清晰度较低,这主要是因为观测区域的AOD整体较高(图4(b))。从MOD04产品中提取该区域范围的AOD,结果表明,该观测区域AOD主要分布在0.3~1.7之间,平均AOD为1.04,区域AOD标准差为0.28,变异系数为26.92%,这说明区域AOD分布离散程度较高,空间差异相对较大。
首先,根据区域AOD平均值,对整副影像进行AOD均值大气校正,通过对比大气校正前(图4(a))和大气校正后的真彩图(图4(c))可以看出,经过大气校正,整幅影像变“亮”且清晰度有所提高,较好地恢复了地表原貌。然后,根据区域AOD分布情况,对整幅影像进行逐像元大气校正,对比真彩图(图4(d))分析可以看出,逐像元大气校正后不仅较大气校正前影像的清晰度大大提高,而且与AOD均值大气校正效果相比,城市、裸地、植被等多种地物的层次表现更加丰富,主要表现为城市、裸地等亮地表亮度相对更高,植被、水体等暗地表亮度相对更低,地物的边缘更加清晰分明,较好地呈现了不同地物的视觉差异特征。
图4 2015年5月25日GF-1卫星WFV1相机大气校正效果对比图
为对大气校正效果进行定量评估,从原始图像、AOD均值大气校正和逐像元大气校正结果中分别选取城镇、裸地和植被3种典型地表区域,为尽量减少数据噪声的影响,每种类型至少选取5万个像元计算每个波段的平均反射率,结果如图5所示。
图5 3种地表类型的不同大气校正结果对比
从图5对比结果可以看出,经过大气校正,城镇、裸地和植被3种典型地物的地表反射率光谱曲线均有所变化。在近红外波段(b4),3种典型地表类型反射率都呈现出逐像元大气校正>AOD均值大气校正>大气校正前;在红波段(b3),城镇和植被反射率呈现大气校正前>AOD均值大气校正>逐像元大气校正结果,裸地反射率则呈现逐像元大气校正>AOD均值大气校正>大气校正前;在绿波段(b2)和蓝波段(b1),城镇和植被反射率呈现大气校正前>AOD均值大气校正>逐像元大气校正结果,裸地反射率则呈现大气校正前>逐像元大气校正结果>AOD均值大气校正。
总体来说,相比较AOD均值大气校正结果,逐像元大气校正一方面更好地消除了城镇和植被大气散射在蓝(b1)、绿(b2)、红(b3)3个可见光波段对地表反射率的增强;另一方面则较好地校正了近红外波段(b4)地表反射的削弱。逐像元大气校正方法在WFV相机4个波段的校正效果明显优于AOD均值大气校正。
表2 典型地物不同大气校正后的NDVI对比
统计3种典型地物的NDVI,发现大气校正后的3种地物NDVI较大气校正前有明显提升(表2),并且通过大气校正进一步提高了不同地物之间NDVI的差异,增强了不同地物之间的辨识度。值得注意的是,逐像元AOD大气校正对地物NDVI的提升效果更为显著,尤其在植被覆盖地区的NDVI增强效果更为突出,这也进一步证明了逐像元大气校正的优势。
根据图4(d)和图5可以看出,一方面在卫星影像的左侧气溶胶较厚的地区经过大气校正后仍存在了部分“雾”状模糊区域,这说明该区域的气溶胶散射影响还未得到完全消除;另一方面大气校正后的地表反射率出现了极少部分的负值,这说明部分像元的气溶胶散射影响有所高估,从而导致校正后的地表反射率偏低。这主要和MOD04产品中AOD的“高值低估、低值高估”的反演误差现象有关,从而导致将MOD04的AOD资料用于大气校正时会发生“高AOD时校正不足、低AOD校正过多”的情况。为进一步分析AOD误差对GF-1卫星WFV相机4个波段大气校正效果的影响,设定表观反射率和观测几何情况(根据图4案例中的GF-1卫星相关参数设定),利用6SV模拟了不同AOD条件下经过大气校正计算的4个波段地表反射率结果,并通过线性拟合分析评估其相对变化率。
图6 GF-1卫星WFV相机反演地表反射率随AOD变化模拟结果
以图4中GF-1卫星WFV1相机观测影像为例,计算整幅影像每个波段的平均表观反射率,然后采用6SV模拟计算在该表观反射率情况下每个波段在不同AOD时大气校正获取平均地表反射率,并采用线性拟合的方法定量估算AOD误差在大气校正过程中给不同波段地表反射率计算结果所带来的误差。从图6可以看出,在一定表观反射率情况下,随着AOD的增高,3个可见光波段大气校正后的地表反射率均呈现逐渐降低的变化特征,而近红外波段地表反射率则呈逐渐升高的变化趋势,这是因为AOD的升高代表大气散射对可见光波段表观反射率的贡献和近红外波段表观反射率的削弱越来越大,大气校正时从表观反射率中扣除大气气溶胶散射的影响也就越来越多,因此,可见光波段和近红外波段的反射率分别呈下降和上升的趋势。根据4个波段的线性拟合结果可以看出,3个可见光波段的线性拟合决定系数较为接近,R2主要在0.86~0.88之间,近红外波段的拟合决定系数明显高于可见光波段,R2在0.97以上,这说明在大气校正时,在近红外波段AOD的变化相比可见光波段更接近线性变化关系。总体上,4个波段的线性拟合决定系数均较高,说明在一定条件下,AOD变化对4个波段大气校正的影响均可用线性关系进行表示,线性拟合的比例系数代表AOD与地表反射率之间的相对变化率。按照MOD04产品AOD约15%的平均误差计算,在一定表观反射率和观测几何条件下,AOD误差给蓝、绿、红和近红外4个波段带来的平均地表反射率误差分别为-0.028、-0.016、-0.007 7和0.009 5。
本文以GF-1卫星WFV1相机观测数据为例,探讨了6S和6SV辐射传输模型对WFV相机4个波段在不同AOD时反应灵敏性,构建了AOD均值大气校正和逐像元大气校正方法,并进行了大气校正误差分析,主要取得以下结论:①通过对6S和6SV模拟GF-1卫星WFV相机探测AOD变化的灵敏性进行比较发现,6SV矢量辐射传输模型模拟效果更好,尤其在高AOD情况下与6S模拟结果产生较大差异,6SV更适用于GF-1卫星WFV相机的准确大气校正;②与传统普遍采用的AOD均值大气校正算法相比,本文提出的逐像元大气校正方法考虑AOD时空变化差异,在还原地物真实地表发射率变化特征方面具有明显优势,不仅能更好地校正大气散射对地表反射率在可见光波段的增强和近红外波段的削弱,还明显提高了不同地表的NDVI差异,能更精准获取GF-1卫星WFV相机的地表反射率;③根据AOD约15%的平均误差计算,在一定表观反射率和观测几何条件下,AOD误差给蓝、绿、红和近红外4个波段带来的平均地表反射率误差分别为-0.028、-0.016、-0.007 7和0.009 5。
本文提出的通过MOD04产品的GF-1卫星WFV相机6SV大气校正方法,考虑了AOD的空间差异,实现了GF-1卫星WFV相机地表反射率更精确的定量反演,该方法在国产的环境一号卫星、资源卫星及GF-1、GF-2、GF-6等卫星影像的大气校正方面具有借鉴意义。为了简化讨论,本文未考虑气溶胶模式、观测几何、地表反射的方向性和邻近效应以及AOD尺度差异等影响,这些简化可能会对大气校正结果带来一定误差。
致谢:本文MOD04资料从NASA戈德太空飞行中心获得,采用的6S和6SV辐射传输模式由Vermote提供,在此一并感谢。