张永红 陈瀚阅 陈宜金 朱利 殷亚秋 杨红艳 侯海倩
(1 中国矿业大学(北京),北京 100083)
(2 福建农林大学,福州 350002)
(3 中国科学院遥感与数字地球研究所,北京 100101)
(4 环境保护部卫星环境应用中心,北京 100094)
(5 中国国土资源航空物探遥感中心,北京 100083)
近年来,核能发电因不会释放温室气体而受到了广泛应用。在核电发电过程中,大量高于自然水温6~11℃的冷却水被排入受纳水体,造成核电基地周围海域水体温度升高,形成了水体热污染[1]。海水表面温度是一项最基本的海洋水文要素,也是海洋—大气系统中一个十分关键的物理量,其变化影响着海洋生物的生存、繁殖、分布特征以及海洋的生态环境等[2],因此对其变化进行监测有着十分重要的意义。传统水表温度获取方法是在区域内布点取样,利用温度测量仪器对采样点进行测量,获得离散点的温度数据,但这种方法实施起来难度大,成本过高且很难实现大面积同步获取数据。热红外遥感技术因其成本低、速度快、资料同步性好、大范围面状观测等优点被广泛应用于海水表面温度的变化监测[3]。2013年2月11日美国宇航局(NASA)发射的Landsat-8卫星携带的热红外传感器(thermal infrared sensor,TIRS)具有两个热红外通道(波段10和波段11),波长范围分别为10.60~11.20μm和11.50~12.50μm,空间分辨率为100m[4],与之前的Landsat系列卫星TM或ETM+传感器相比,热红外波段由单通道变为双通道;与NOAA系列卫星的AVHRR数据和MODIS数据相比,空间分辨率显著提高;与国内环境与灾害监测卫星HJ-1B热红外数据相比,波段数量增多且空间分辨率较高,因此,Landsat-8数据源具备较高的优势,可作为国内卫星监测核电基地周围海表温度的有效补充和验证,有必要对其进行进一步应用研究。
Landsat-8/TIRS有两个热红外通道,海表温度遥感反演不仅可以利用单通道成像数据,还可以利用双通道成像数据。基于单通道的反演算法有辐射传输方程法(radiative transfer model,RTM)、覃志豪等提出的单窗算法(mono-window model,MW)、Jimenez-Munoz和Sobrino提出的普适性单通道算法(JM&S)[5-7]。基于双通道的有Price最早针对NOAA卫星的第4、5通道提出的Price劈窗算法,Becker-Li在Price劈窗算法基础上改进的Becker-Li劈窗算法,Wan和Dozier将Becker-Li劈窗算法引用到MODIS传感器并进行改正的Wan劈窗算法,覃志豪在Wan劈窗算法基础上针对NOAA/AVHRR数据改进的Qin劈窗算法(Qin split window,QSW)[8]。目前,基于Landsat-8/TIRS数据,对核电基地周边海域海表温度遥感反演算法适用性比对分析的研究较少,因此本文选择辽宁红沿河核电基地周围海域为研究区,采用辐射传输方程法、单窗算法、普适性单通道法和Qin劈窗算法分别进行算法参数区域适用性校正和海表温度反演,并结合海面实测数据对其反演精度进行比对分析,为后续Landsat-8/TIRS数据应用于红沿河核电基地遥感监测和我国其他核电基地遥感监测提供参考和借鉴。
红沿河核电基地海域位于渤海辽东湾东海岸大连市瓦房店温坨子村。该区属于温带季风气候区,以大陆性气候为主,季风显著,降水集中,四季分明,多年平均年降水量为540.3mm,大气稳定度以中性和稳定为主,大气扩散条件良好,海域生物种类繁多。其中斑海豹是国家二级保护动物,是鳍脚类动物在中国繁殖的唯一种类,而斑海豹自然保护区则位于核电厂附近海域,斑海豹每年秋冬季节洄游到该海域,登上浮冰繁殖,春季冰融后分散在沿岸觅食,在此期间如果浮冰融化,将不利于斑海豹生长繁殖。
影像数据选取 2014年 6月 11日 Landsat-8卫星于红沿河核电基地海域过境影像,图幅号为LC81200322014162LGN00。对该影像首先进行波段合并、裁剪以及分别对波段10和波段11辐射定标,获取星上辐亮度、计算亮温等数据预处理;再根据4种温度反演算法获取海表温度,最后与实测数据进行比对分析。实测数据采用温盐深仪CTD于卫星过境前后30min内测量,表面测量深度在0.5m左右。
本文选择了辐射传输方程法、单窗算法、普适性单通道法和Qin劈窗算法4种温度反演算法进行比对分析,由于前三种算法基于一个热红外通道,波段10 和波段11都可以,但是波段10的光谱范围处于大气窗口内大气吸收较低区域,大气透过率高,反演效果较好[9],因此选用波段10作为前三种算法反演的通道选择。
(1)辐射定标
根据Landsat-8/TIRS热红外数据头文件提供的信息,利用绝对定标系数将灰度值(DN)图像转换为辐亮度图像,星上辐亮度Lλ计算公式为
式中 DN为灰度值;h为偏移量,取值0.1;Gain为绝对定标系数增益,取值3.342 0×10–4。
(2)亮温计算
在获取星上辐亮度的基础上,利用Planck函数反函数近似式计算亮度温度,计算公式为
式中 Ti为波段10或波段11的像元亮度温度(热力学温度);k1和k2为发射前预设的常量,通过头文件查询(波段10,k1=774.89,k2=1 321.08;波段11,k1=480.89,k2=1 201.14)。
2.2.1 辐射传输方程法
根据热红外辐射传输模型推导,在晴空条件下不考虑多次散射,假设水平均匀大气,对于一定海表温度Ts,海表的黑体辐射 B(λ, Ts)为
式中 Lsensor为星上辐亮度,由卫星影像定标获取;Lup和Ldown分别为大气上下行辐射,τ为大气透过率,这三个参数可在NASA官网中输入相关参数后获取;ε为地表发射率(对波段10的水体取值为0.992)[10]。
2.2.2 单窗算法
覃志豪等通过大气平均作用温度Ta将大气辐射简化,根据热辐射传输方程以及一系列假设,建立了一个适用于Landsat/TM6的地表温度计算算法,如式(4)所示:
式中 C、D是中间变量,其中,C=τε, D=(1–τ)[1+(1–ε)τ];a、b为回归系数;Tsensor为星上亮温。
(1)大气平均作用温度Ta估算
不同层的大气温度往往是不可知的,覃志豪等人提出了在四种标准大气廓线的模式下根据近地表空气温度To来估计大气平均作用温度Ta,如表1所示。
表1 大气平均作用温度计算模型Tab. 1 The calculation model of average atmospheric temperature
(2)大气透过率τ估算
大气透过率τ可以根据大气总水汽含量w来估算,通常做法是基于TIGR(thermodynamic initial guess retrieval)大气廓线数据,结合MODTRAN软件模拟τ与w之间的关系。本文基于中纬度夏季大气模式,以w分段线性拟合,建立针对Landsat-8波段10的大气透过率估算方程:
大气总水汽含量w可根据美国国家环境预报中心(NCEP)探空数据计算。本文从NASA官网获取卫星过境时研究区 NCEP探空数据[11],以此获取上空相对湿度、气压、气温等资料,并通过式(6)计算大气总水汽含量w。
式中 ρ为液态水密度;gn为标准重力加速度;T0为热力学标准态温度,T0=273.135K;p0为标准大气压,p0=1 013hPa;T、p和q为规定层大气的热力学温度、气压和上空相对湿度;pz为海拔高度z处的大气压。
(3)回归系数估算
单窗算法中对普朗克公式按照泰勒级数线性展开,定义了一个温度参数L。L的数值与温度T有很好的线性关系,如式(7)所示。
针对Landsat-8波段10光谱响应函数,在0~70℃温度变化范围内线性拟合L和T(如式(8)所示),计算得出a=–66.323,b=0.446 4,相关系数R2=0.999 4。
2.2.3 普适性单通道算法
Jiménez-Muñoz和Sobrino’s对普朗克函数在某个温度值附近作一阶泰勒展开,并以高斯三角滤波函数作为热红外波段通道响应函数对大气效应进行模拟,得到地表温度Ts反演算法,
式中 ψ1、ψ2、ψ3为大气影响因子,是大气总水汽含量w的函数;γ,δ为中间变量,计算公式为
其中 c1=1.191 04×108W·m–2·sr–1·μm–1,c2=1.438 8×104μm·K;λeff是中心波长,由波长λ与通道响应函数f(λ)卷积而得。
在标准大气下,Landsat-8波段10的通道范围较窄,选用通道中心波长10.9μm作为等效波长,利用MODTRAN模拟结果拟合得到大气影响因子ψ1、ψ2、ψ3关于大气总水汽含量w的函数关系:
2.2.4 Qin劈窗算法
覃志豪通过大气向下热辐射的近似解和Planck辐射函数的线性化对前人的劈窗算法进行推导,得到适用于NOAA卫星AVHRR探测器的两个热红外通道的温度反演算法(Qin劈窗)[12]。该算法仅需两个参数(大气透过率和地表发射率)对地表温度进行反演,本文根据Landsat-8卫星第10、11通道与AVHRR卫星第4、5通道的热辐射相似性,对Qin劈窗算法(QSW)进行修订[14],如式(11)所示。
式中 Ts为海表温度;T10和T11分别为波段10和波段11的亮度温度;A0、A1、A2分别为中间变量,计算公式分别为:
其中 Ci,Di为波段10或11对应的中间变量,, Di=[1 - τi( θ) ][1+(1 - εi)τi(θ)];εi为波段10(或波段11)对应的地表发辐射率(波段11的水体取值为0.989);大气透过率τi(θ)可由单窗算法中求取大气透过率方法获取;系数a10、b10、a11、b11可通过单窗算法中求取参数a、b方法获取。
根据修订的模型参数,四种反演算法得到的海表温度反演结果如图1所示,灰度颜色由深到浅表示温度变化由高到低。从图1中可以看出,从海洋到陆地温度逐渐升高,QSW法整体呈现浅色调,反演温度最高;RTM法较之加深;MW算法呈现深色调,反演温度最低,JM&S法与MW算法相近。
利用海表温度实测结果减去四种算法的反演结果,绘制曲线如图2所示。从图中可以看出,四种算法反演出的结果总体变化趋势比较接近,其中RTM法所得结果与实测值最为相近,绝对误差较小;MW法和JM&S法反演结果相近,整体较实测值偏小,绝对误差偏大;QSW法反演结果较实测值整体偏大,绝对误差也较大。
从定量方面对各算法误差进行分析统计,如表2所示。可以看出,RTM算法平均偏差接近0,为0.27℃。MW算法平均偏差正向最大,为3.58℃,QSW算法平均偏差负向最大,为–2.86℃。从均方根误差来看,RTM算法最小,为0.61℃,QSW算法次之,为2.93℃,JM&W算法较高,为3.26℃,MW算法最大,为3.64℃。
综上可知,四种算法中RTM反演的海表温度精度最高,QSW算法次之,MW算法最差。
图1 研究区海表温度反演结果Fig. 1 The results of temperature retrieval
图2 不同算法绝对误差比对Fig. 2 The absolute error graph of four kinds of algorithm results
表2 不同反演算法误差分布Tab. 2 The error distribution of four kinds of retrieval algorithm℃
图3反映了四种算法反演海表温度与实测海表温度的关系。对比各算法反演值与实测值的线性关系发现,RTM法反演值与实测值一致性最好,反演过程中,在 NASA官网输入卫星过境时研究区位置信息、时间信息以及气象信息等,由其计算Landsat-8波段10的大气上行辐射、大气下行辐射、大气透过率等参数,这些参数利用准实时大气廓线进行模拟,精度高。MW法反演值小于实测值,在反演过程中使用标准大气廓线进行大气平均作用温度经验公式推导,增加了该算法的误差[14];同样,JM&S法反演值小于实测值,在反演过程中,精确的大气水汽含量是反演的关键,文献[15]研究表明若大气水汽含量较大,反演误差也会随之增大。
图3 各算法反演SST与实测SST线性关系Fig. 3 The linear relationship of sea surface temperature between inversion values and actually measured values
QSW法反演值大于实测值,根据文献[16-17],劈窗算法是较单通道算法发展较为成熟、精度较高的温度反演算法。由于Landsat-8/TIRS传感器接收到了除正常电磁波信号以外的其他散射光信号,该假信号造成了波段10和波段11影像产生条带和辐射校正误差(波段11的影响程度是波段10的两倍多)[18]。NASA地面校正团队通过扫描亮度较高的月球,以此寻找散射光来源,除去了条带影响,减少了辐射校正误差,使得波段10辐射误差在1K要求范围之内,但波段11仍处于1.75~4.4K的较高误差范围[19],其不确定性增大了劈窗算法反演海表温度的误差。
基于2014年6月11日红沿河核电基地地区的Landsat-8/TIRS数据,分别对RTM法、MW算法、JM&S法、QSW算法等进行海表温度反演参数修订与数据处理。反演结果与星地同步实测温度对比分析表明:RTM法反演温度与实测温度相近,MW算法和JM&S法低于实测温度,QSW算法高于实测温度;另一方面,从均方根误差统计来讲,RTM法误差最小,QSW法次之,MW算法和JM&S法误差最大。因此,基于Landsat-8/TIRS热红外数据进行红沿河核电基地遥感监测,辐射传输方程法温度反演精度最高,监测效果最好,可业务化应用于红沿河核电基地周边海域海表温度遥感反演,为核电温排水排放监测提供技术保障。
(
)
[1] 刘永叶, 陈晓秋. 核电厂温排水热影响研究的建议[J]. 辐射防护通讯, 2011, 31(6): 20-23. LIU Yongye, CHEN Xiaoqiu. Suggestion of the Study on Thermal Impact of Thermal Discharge from NPPs[J]. Radiation Protection Bulletin, 2011, 31(6): 20-23. (in Chinese)
[2] 张彩, 周世健, 朱利. 基于环境星HJ-1B数据反演大亚湾海表温度的三种算法[J]. 东华理工大学学报, 2013, 36(1): 88-92. ZHANG Cai, ZHOU Shijian, ZHU Li. Three Algorithms of Sea Surface Temperature Inversion of Daya Bay Based on Environmental Satellite HJ-1B Data[J]. Journal of East China Institute of Technology, 2013, 36(1): 88-92. (in Chinese)
[3] 朱利, 赵利民, 王桥, 等. 核电站温排水分布卫星遥感监测及验证[J]. 光谱学与光谱分析, 2014, 34(11): 3079-3084. ZHU Li, ZHAO Limin, WANG Qiao, et al. Monitoring the Thermal Plume from Coastal Nuclear Power Plant Using Satellite Remote Sensing Data: Modeling and Validation[J]. Spectroscopy and Spectral Analysis, 2014, 34(11): 3079-3084. (in Chinese)
[4] 张玉君. Landsat-8简介[J]. 国土资源遥感, 2013, 25(1): 176-177. ZHANG Yujun. Landsat-8 Brief Introduction[J]. Remote Sensing for Land & Resources, 2013, 25(1): 176-177. (in Chinese)
[5] 韩启金, 傅俏燕, 潘志强, 等. 利用HJ-1B星热红外遥感图像研究城市热岛效应[J]. 航天返回与遥感, 2012, 33(1): 67-74. HAN Qinjin, FU Qiaoyan, PAN Zhiqiang, et al. Study on Urban Heat Island Effect with HJ-1B Thermal Infrared Remote Sensing Image[J]. Spacecraft Recovery & Remote Sensing, 2012, 33(1): 67-74. (in Chinese)
[6] 覃志豪, Zhang M H, Karnieli A, 等. 用陆地卫星TM6数据演算地表温度的单窗算法[J]. 地理学报, 2001, 56(4): 456-466. QIN Z H, Zhang M H, Karnieli A. Mono-window Algorithm for Retrieving Land Surface Temperature from Landsat TM6 Data[J]. Acta Geographica Sinica-Chinese Edition-, 2001, 56(4): 466-474.
[7] Jimenez-Munoz J C, Sobrino J A. A Generalized Single-channel Method for Retrieving Land Surface Temperature from Remote Sensing Data[J]. Journal of Geophysical Research: Atmospheres (1984-2012), 2003, 108(D22).
[8] Qin Z, Karnieli A. Progress in the Remote Sensing of Land Surface Temperature and Ground Emissivity Using NOAA-AVHRR Data[J]. International Journal of Remote Sensing, 1999, 20(12): 2367-2393.
[9] Jimenez-Munoz J C, Sobrino J, Skokovic D, et al. Land Surface Temperature Retrieval Methods from Landsat-8 Thermal Infrared Sensor Data[J]. Geoscience and Remote Sensing Letters, IEEE, 2014, 11(10): 1840-1843.
[10] 杨槐. 从Landsat-8影像反演地表温度的劈窗算法研究[J]. 测绘地理信息, 2014, 39(4): 73-77. YANG Huai. Research of Split Window Algorithm for Retrieval of Land Surface Temperature from Landsat-8[J]. Journal of Geomatics, 2014, 39(4): 73-77. (in Chinese)
[11] USGS. Atmospheric Correction Parameter Calculator[EB/OL].[2014-02-18]. http://atmcorr.gsfc.nasa.gov/.
[12] Qin Z, Dall’Olmo G, Karnieli A, et al. Derivation of Split Window Algorithm and Its Sensitivity Analysis for Retrieving Land Surface Temperature from NOAA-advanced very High Resolution Radiometer Data[J]. Journal of Geophysical Research: Atmospheres (1984-2012), 2001, 106(D19): 22655-22670.
[13] Yang L, Cao Y G, Zhu X H, et al. Land Surface Temperature Retrieval for Arid Regions Based on Landsat-8 TIRS Data: a Case Study in Shihezi, Northwest China[J]. Journal of Arid Land, 2014, 6(6): 704-716.
[14] 白洁, 刘绍民, 扈光. 针对TM/ETM+遥感数据的地表温度反演与验证[J]. 农业工程学报, 2008, 24(9): 148-154. BAI Jie, LIU Shaomin, HU Guang. Inversion and Verification of Land Surface Temperature with Remote Sensing TM/ETM+ Data[J]. Transactions of the CSAE, 2008, 24(9): 148-154. (in Chinese)
[15] Jimenez-Munoz J C, Cristobal J, Sobrino J A, et al. Revision of the Single-channel Algorithm for Land Surface Temperature Retrieval from Landsat Thermal-infrared Data[J]. Geoscience and Remote Sensing, IEEE Transactions on, 2009, 47(1): 339-349.
[16] 毛克彪, 覃志豪, 施建成. 用 MODIS影像和劈窗算法反演山东半岛的地表温度[J]. 中国矿业大学学报, 2005, 34(1): 46-50. MAO Kebiao, QIN Zhihao, SHI Jiancheng. Retrieval of Land Surface Temperature of Shandong Peninsula Using MODIS Image and Split Window Algorithm[J]. Journal of China University of Mining& Technology, 2005, 34(1): 46-50. (in Chinese)
[17] 朱利, 顾行发, 陈良富, 等. 高分辨率红外相机单窗与劈窗陆表温度反演精度分析研究[J]. 红外与毫米波学报, 2008, 27(5): 346-353. ZHU Li, GU Xingfa, CHEN Liangfu, et al. Comparison of LST Retrieval Precision Between Sigle-channel and Split-Window for High-resolution Infrared Camera[J]. Journal of Infrared and Millimeter Waves, 2008, 27(5): 346-353. (in Chinese)
[18] Montanaro M, Gerace A, Lunsford A, et al. Stray Light Artifacts in Imagery from the Landsat-8 Thermal Infrared Sensor[J]. Remote Sensing, 2014, 6(11): 10435-10456.
[19] USGS. Landsat Headlines[EB/OL]. [2013-11-14]. http://landsat.usgs.gov/calibration_notices.php.