0 引言

自1988年Gomy等发现1984年加兹利地震震前卫星热红外异常后（张帆等，2016），有关卫星热红外的研究逐渐开展起来，研究结果已在震后紧急救援、灾情监测、灾后重建等方面得到应用（叶秀薇等，2010）。

20世纪80年代末，马瑾在国内最先开展了卫星热红外与地震间相关性的研究（张治广等，2017）；在此基础上，强祖基与徐秀登等合作开展了热红外遥感研究（徐好民，2002）。随着遥感技术的发展，卫星数据质量不断提高，利用气象卫星资料，获取海面温度的技术趋于成熟，海面温度遥感反演技术的成功使得地表温度遥感反演成为遥感领域又一研究热点（黄启厅，2010）。近年来，地表温度异常、植被覆盖度变化与地震异常间相关性的研究逐渐兴起。胡光庭等（2017）分析了阳江市地表温度变化特征，认为植被覆盖度的增加能很好地降低地表温度；王蕾（2014）等基于遥感影像发现，植被覆盖度对缓解城市热岛效应具有一定作用，说明地表温度与植被覆盖度间具有较好的的负相关性；庞光辉等（2016）对沈阳市植被覆盖度降温效应进行了分析，发现植被覆盖度每增加10%，地表温度可降低0.71℃。前人研究结果大多适用于城市生态环境规划和建设，且反演方法较简单，对地表物未进行详细分类而直接提取地表温度数据，且未统计地表温度平均值及标准差。本文采用辐射传输方程算法对研究区域地表温度进行了反演，为提高反演结果的准确性，使用最大似然分类算法对地表物进行了耕地、林地、草地、水体、建设用地、未利用土地的详细分类，依据分类结果及反演结果，分别计算了6种地表物所在区域内的地表温度、植被覆盖度的平均值、标准差，分析总结了研究区域地表温度与植被覆盖度间的相关性。

1 研究区及数据

乌兰浩特市地理坐标为45°22′—46°18′N、121°51′—122°20′E，总面积2 353.5 km²，属中温带半干旱季风气候，四季分明，气温的年差、日差都较大，年平均气温5℃，年降水量为440 mm左右。乌兰浩特地震台位于内蒙古自治区兴安盟乌兰浩特，是兴安盟唯一有人值守的测震数字化台站，承担着兴安盟区域内地震监测预测任务（包宝小等，2019）。本文选取过境时间为2010、2015、2020年3年的9月11日2时乌兰浩特地震台监测区域的Landsat 5、Landsat 8遥感影像数据（云量较少，符合使用需求），进行辐射定标、大气校正预处理，使用乌兰浩特矢量图进行矢量裁剪。利用影像数据第1、2、3、4、5、7波段提取NDVI（归一化植被指数）及土地利用分类，利用Landsat 5的第6波段及Landsat 8的第10波段反演地表温度数据。

2 原理 2.1 辐射定标

辐射定标是将传感器记录的无量纲DN值转换成具有实际物理意义的大气顶层辐射亮度或反射率。其原理是建立数字量化值与对应视场中辐射亮度之间的定量关系，以消除传感器本身产生的误差（央金卓嘎等，2019）。

2.2 辐射传输方程法

卫星传感器接收到的热红外辐射亮度L_λ由3部分组成：大气向上辐射亮度L↑、地面的真实辐射亮度经过大气层之后到达卫星传感器的能量、大气向下辐射到达地面后反射的能量L↓。卫星传感器接收到的热红外辐射亮度L_λ的表达式可写为（辐射传输方程）

$ L_\lambda=\left[\varepsilon B\left(T_{\mathrm{s}}\right)+(1-\varepsilon) L \downarrow\right] \tau+L \uparrow $

(1)

式中，ε为地表比辐射率；T_s为地表真实温度；B(T_s)为通过反映物体温度和辐射强度之间关系的普特朗克公式推算所得的温度T相对应的黑体辐射亮度；τ为大气在热红外波段的透射率；L↑、L↓为大气的向上、向下辐射亮度（崔耀辉，2018）。由此可得B(T_s)的计算公式如下

$ B\left(T_{\mathrm{s}}\right)=\left[L_\lambda-L \uparrow-\tau(1-\varepsilon) L \downarrow\right] / \tau \varepsilon $

(2)

其中，L↑、L↓可通过影像的中心纬度与成像时间在网站（http://atmcorr.gsfc.nasa.gov）获取（张菊等，2020）。

运用普朗克公式的反函数，可推导出亮度温度T_R

$ T_R=\frac{K_2}{\ln \left[\frac{K_1}{B\left(T_{\mathrm{s}}\right)}+1\right]} $

(3)

其中，K₁、K₂的值可在影像MLT.txt文件中获得。

2.3 地表比辐射率及植被覆盖度计算

地表比辐射率又称发射率，是指在同一温度下地表发射的辐射量与一黑体发射的辐射量间的比值，是利用热红外遥感获取地表温度的必不可少的参数。针对Landsat影像的第10段热红外波段，使用Sobrino提出的NDVI阈值法计算地表比辐射率（齐广慧等，2020）

$ \varepsilon=0.004 P_{\mathrm{v}}+0.986 $

(4)

其中，P_v为植被覆盖度，其公计算式如下

$ P_{\mathrm{v}}=\left(\mathrm{NDVI}-\mathrm{NDVI}_{\mathrm{s}}\right) /\left(\mathrm{NDVI}_{\mathrm{v}}-\mathrm{NDVI}_{\mathrm{s}}\right) $

(5)

其中，NDVI为归一化植被指数；NDVI_s为完全是裸土或无植被覆盖区域的NDVI；NDVI_v为完全植被覆盖的像元的NDVI。一般取经验值，NDVI_s =0.05，NDVI_v=0.70（蔡贤等，2020）。

2.4 最大似然分类

为提高地表温度反演的准确性，使用最大似然分类方法对影像进行地表物的分类。在2类或多类判别中，假定各类分布函数为正态分布，并选择训练区，计算各待分类样区的归属概率而进行分类的图像分类方法称为贝叶斯分类方法，本文根据贝叶斯准则对遥感影像进行分类。贝叶斯定理：假设B₁、B₂，…互斥且构成一个完全时间，A伴随它们出现，已知它们分别发生的先验概率P(B_i)（i = 1，2，…）及A的条件概率P(B_i)，则可以得到事件A的后验概率P(B_i|A)。由概率乘法公式

$ P(A B)=P(A) * P(B \mid A)=P(B) * P A(A \mid B) $

(6)

可导出贝叶斯定理公式

$ P(A \mid B)=\frac{P(B \mid A) * P(A)}{P(B)} $

(7)

设有s个类别，用ω₁、ω₂、…、ω_s来表示，每个类别发生的概率（先验概率）分别为P(ω₁)、P(ω₂)、…、P(ω_s)。设有未知类别的样本X，其类条件概率分别为

$ P\left(X \mid \omega_1\right), P\left(X \mid \omega_2\right), \cdots, P\left(\mathrm{X} \mid \omega_s\right) $

(8)

则根据贝叶斯定理可以得到样本X出现的后验概率为

$ P\left(\omega_i \mid X\right)=\frac{P\left(X \mid \omega_i\right) * P\left(\omega_i\right)}{P(X)}=\frac{P\left(X \mid \omega_i\right) * P\left(\omega_i\right)}{\sum_{i=1}^s P\left(X \mid \omega_i\right) * P\left(\omega_i\right)} $

(9)

此时，以样本X出现的后验概率作为判别函数来确定样本X的所属类别，其分类准则为：如果P(ω_i|X) = max_{1 < j < s}P(ω_j|X)，则X∈ω_i；通过观测样本X，把先验概率P(ω_i)转化为后验概率P(ω_i|X)，并以后验概率为最大的原则确定样本X的所属类别。

3 计算结果与分析

利用影像数据的红光波段及近红外波段提取了2010、2015、2020年的归一化植被指数（NDVI）（图 1）。NDVI的取值范围-1—1，当NDVI值＜0时，表示水域，即乌兰浩特及科尔沁右翼前旗区域间的归流河；当NDVI值为0时，表示地表被裸土及岩石覆盖，为建筑用地，即乌兰浩特主城区区域；当NDVI值为1时，表示植被覆盖丰富，即草地、林地及耕地、主城区四周区域。由图 1可见，2010—2020年NDVI低值区域没有太大变化，高值区域明显向四周扩散。

为反演结果的准确性，本文使用监督分类方法中的最大似然分类方法将影像进行6种地表物分类，分为耕地、林地、草地、水体、建设用地及未利用地，分类结果如图 2所示。

遥感影像Landsat TM/ETM+的第6、10、11波段属于热红外波段，可用于地表温度反演，地表温度反演算法主要有辐射传输方程法、单窗算法、单通道算法、基于影像的反演算法（张菊等，2020）。本文使用辐射传输方程法对研究区域影像进行地表温度反演，反演结果如图 3所示。

依据分类结果及反演结果影像，分别统计了6种地表物所在区域内的地表温度及植被覆盖度的平均值与标准差，统计结果如表 1、2所示。

由表 1、2可见，水体因自身的物理特征，温度最低；其次是耕地，温度均值为29.5℃，标准差也相对较低；再次为耕地及林地；草地地表温度均值最高，说明不同地表物之间地表温度差异较大。建设用地与未利用土地相比，未利用土地地表温度较高，建设用地上虽然生长的植物较少，但乌兰浩特市绿化较好，所以建设用地地表温度比未利用土地低。因建设用地的绿化结构较复杂，由表 1、表 2可见，建设用地地表温度及植被覆盖度的标准差较大。地表温度较低的地表物类型，其植被覆盖度相对较高；而地表温度较高的地表物类型，植被覆盖度则较低，这说明植被有很好的降温作用。所以，进行地表温度映震异常分析时要多方面考虑植被覆盖度的影响。

4 结论与讨论

选取乌兰浩特地震台监测区域2010、2015、2020年3年的9月11日2时的Landsat 5、Landsat 8遥感影像数据，经过辐射定标、大气校正、裁剪、拼接等预处理，提取NDVI（归一化植被指数），计算植被覆盖度，使用最大似然分类法进行土地利用分类，使用Landsat 5的第6波段及Landsat 8的第10波段反演了地表温度。研究得到如下结论：①将林地、草地、耕地分为一类，按照地表温度进行排序，从低到高为耕地、林地、草地，建设用地及未利用土地分为一类，未利用土地地表温度高于建设用地；②按照植被覆盖度从低到高排序结果为草地、林地、耕地；③因水体自身特征，地表温度及植被覆盖度都很低，除水体之外，建设用地地表温度及植被覆盖度标准差高于其他地表物；④由不同地表物反演的地表温度有明显差异，不同地表物植被覆盖度亦各异；⑤地表温度较低的地表物类型植被覆盖度相对较高，地表温度较高的相对较低，说明对于不同地表物类型植被均有降温作用。

综上所述，地表温度异常不一定为地震短临异常，需综合考虑如下干扰因素：当地表温度出现增温现象时，首先应考虑是否植被覆盖度降低所致，如表 1中2015年与2020年建设用地地表温度均值出现上升现象，但并没出现地震事件，而是市政规划将草地变为建设用地所致。当地表温度出现降温现象时，需实地取证，以明确是否增添植被而导致地表温度下降的，如表 1中草地、林地、耕地地表温度逐年下降，但并没出现地震事件，而是因为植被覆盖度增加了。当地表温度斑块性异常时，需根据市政规划文件，实地考察是否建设用地进行了重新规划改造所致。本文给出的2010、2015、2020年3年的9月地表温度、植被覆盖度的均值及标准差，所得结果可为研究地震热红外异常工作提供有效的数据支持及科学依据。