本发明涉及遥感技术领域，尤其涉及的是一种山地地表温度遥感反演方法。
背景技术：
地表温度是描述地表能量分布情况以及地表与大气相互作用的关键因素，也是研究地表物质能量动态平衡的重要参数，被广泛应用于灾害评估、蒸散发估算、气候研究和土壤水分监测等方面。目前的地表温度反演方法大多针对平坦地面，开展复杂地形区的研究甚少。我国五大典型地形包括平原、山地、高原、盆地和丘陵，其中平原面积仅占12％。在复杂地形区，由于地形的起伏造成地表结构发生变化，这种变化引起地面与传感器的几何关系发生改变，进而导致传感器接收到的地面辐射能量的组成部分发生变化。现如今，越来越多高空间分辨率的卫星数据盛行，地形影响在高空间分辨率数据中尤为突出，这使得我们不得不考虑地形起伏导致的辐射能量组成部分的变化。
现有地表温度反演方法主要适用于平坦地表，针对复杂地形区的研究甚少。然而，山地地形的复杂几何结构使得地表与大气之间的热辐射传输方式发生变化，这使得基于常规热辐射传输方程的地表温度反演方法不能适用于山地地表温度反演。
因此，现有技术存在缺陷，需要改进。
技术实现要素：
本发明所要解决的技术问题是针对现有技术的不足提供一种山地地表温度遥感反演方法。
本发明的技术方案如下：
一种山地地表温度遥感反演方法，包括以下步骤：
步骤1：下载搭载在landsat8卫星上的tirs传感器的星上辐亮度产品，以及oli传感器的地表反射率产品，并进行数据预处理；
步骤2：地表发射率计算；
步骤3：地形参数计算；
下载地表高程dem数据，地形参数计算过程包括：
(1)计算坡度和坡向
坡度和坡向采用三阶反距离平方权差分法计算：
式中，s和a分别为坡度和坡向，slopens为以坡面目标点高程为中心的3×3个像元内的南北走向趋势，slopewe为东西走向的趋势，分别通过公式(12)和(13)计算：
slopewe＝(hj 1,i-1 2hj 1,i hj 1,i 1-hj-1,i-1-2hj-1,i-hj-1,i 1)/(8×cellsize)(22)
slopens＝(hj-1,i 1 2hj,i 1 hj 1,i 1-hj-1,i-1-2hj,i-1-hj 1,i-1)/(8×cellsize)(23)
式中，h为像元的地表高程，i和j分别为像元在地表高程图像中的行号和列号，cellsize为像元的空间分辨率；
(2)计算天空可视因子
天空可视因子表示山地地表一个点在半球空间内对天空的可视程度，即该点接收到的大气下行辐射与其在平坦地表未被遮挡时接收到的来自半球空间内的大气下行辐射的比值，取值范围为0到1，可以通过公式(14)计算：
式中，s为坡度，a为坡向，φ为方位角，hφ为水平张角，它表示当前目标像元所在水平面的法线与当前目标像元和周围邻近像元之间连线的夹角；将方位角离散为均匀分布的16个方向计算天空可视因子，这16个方向分别为0°,30°,45°,60°,90°,120°,135°,150°,180°,210°,225°,240°,270°,300°,315°,和330°；
步骤4：大气参数计算；大气参数包括大气透过率、大气上行辐射和大气下行辐射；
步骤5：山地地表温度反演；
在热红外波谱范围内，常规的热辐射传输方程表示为：
ltoa＝εb(ts)τ (1-ε)ldownτ lup(25)
式中，ltoa为大气层顶接收到的辐亮度，ε为地表发射率，b为普朗克函数，ts为地表温度，τ为大气透过率，ldown为大气下行辐射，lup为大气上行辐射；
山地热辐射传输方程表示为：
ltoa＝εb(ts)τ (1-ε)ldownτvd (1-ε)ladjτ lup(26)
式中，vd为天空可视因子，ladj为邻近地形辐射，表示为：
式中，下标a和b分别代表目标像元和邻近像元，lb为邻近像元的辐亮度，ta和tb分别为目标像元和邻近像元的法线与目标像元和邻近像元的中心点连线的夹角，dsa和dsb分别为目标像元和邻近像元的坡面面积，rab为目标像元和邻近像元之间的距离，n为邻近像元b的总个数；
基于山地热辐射传输方程反演山地地表温度表示为：
式中，tsm为反演的山地地表温度，b-1为普朗克函数的逆运算；
基于常规的热辐射传输方程反演地表温度表示为：
式中，tsf为反演的平坦地表温度；
具体的山地地表温度反演过程包括：
(1)使用步骤1中的landsat8第10通道的辐亮度、步骤2中计算的地表发射率，步骤4中计算的大气参数，通过公式(19)逐像元反演得到平坦地表温度tsf，并将此时的地表温度作为邻近像元地表温度的初始值；
(2)将邻近像元地表温度的初始值输入公式(17)得到邻近像元辐亮度，并将步骤2中计算的地表发射率，步骤3中计算的坡度、坡向和天空可视因子，步骤4中计算的大气参数，以及步骤1中的landsat8第10通道的辐亮度，同时输入公式(18)，逐像元反演得到山地地表温度tsm，并将此时的地表温度作为邻近像元的地表温度，用于迭代计算邻近像元辐亮度；
(3)重复(2)，直到前后两次邻近像元地表温度的差值小于0.01k或者迭代次数大于3时输出结果，得到最终的山地地表温度。
所述的山地地表温度遥感反演方法，所述步骤1中数据预处理包括:
(1)将星上辐亮度产品中的计数值转换为星上辐亮度：
l＝gain×dn offset(30)
式中，l为星上辐亮度；dn为计数值；gain和offset为增益和偏移；landsat8第10通道的增益和偏移分别为0.00033420和0.1；
(2)将地表反射率产品中的计数值转换为地表反射率：
ρ＝a×dn b(31)
式中，ρ为地表反射率，dn为计数值，a和b为转换系数；landsat8第4通道和第5通道的系数a和b分别为0.0000275和-0.2；
(3)计算归一化植被指数ndvi：
式中，ρred和ρnir分别为红光波段和近红外波段的地表反射率；
(4)计算植被覆盖度：
式中，pv为植被覆盖度，ndvi为当前像元的归一化植被指数，ndvimin和ndvimax分别为ndvi最小值和最大值，ndvimin取值为0.05，ndvimax取值为0.85。
所述的山地地表温度遥感反演方法，所述步骤2中，利用asterged产品计算landsat8第10通道的地表发射率，该产品是利用2000年到2008年所有的晴空aster数据通过温度与发射率分离算法反演的地表温度和发射率平均计算得到，包括aster通道10-14五个热红外通道发射率的平均值和标准偏差，以及其他数据集；landsat8地表发射率计算过程包括：
(1)从asterged产品中提取出ndvi、第13和14通道的dn值；根据公式(5)将dn值转换为asterged的地表发射率：
ε＝c×dn d(34)
式中，ε为asterged的地表发射率，dn为计数值，c和d为转换系数；aster第13和14通道的转换系数c和d分别为0.001和0；
(2)假定像元的地表发射率可以简单地表示为植被与裸土比辐射率的加权值：
εi＝εvipv εsi(1-pv)(35)
式中，εi为通道i的地表发射率，εvi为通道i的植被组分发射率，εsi为通道i的裸土组分发射率，pv为植被覆盖度；
(3)基于asterged的地表发射率和公式(6)，计算aster第13和14通道的裸土组分发射率：
式中，上标a代表aster传感器，i代表aster第13和14通道，
(4)从aster地物波谱库中挑选41条土壤反射率曲线，分别利用aster第13和14通道以及landsat8第10通道的光谱响应函数与反射率波谱曲线进行卷积计算，得到每个通道卷积后的土壤发射率；以aster第13和14通道卷积后的土壤发射率作为自变量，以landsat8第10通道卷积后的土壤发射率作为因变量，拟合得到三个通道土壤发射率之间的回归方程：
式中，
(5)landsat8地表发射率通过土壤组分发射率和植被组分发射率加权计得到：
式中，
所述的山地地表温度遥感反演方法，所述步骤4中，下载ncep大气廓线数据；大气参数计算过程包括：
(1)根据landsat8卫星过境的utc时间，选取最近两个时刻的大气廓线数据；
(2)提取ncep大气廓线31个气压层的高度、温度、相对湿度数据；
(3)将大气廓线输入大气辐射传输模型modtran，计算每条ncep大气廓线对应的大气参数，包括大气透过率、大气上行辐射和大气下行辐射；
(4)根据像元的地表高程，将距离该像元最近的ncep四个网格点的大气参数进行线性插值，插值到同一高程上；
(5)根据像元中心的经纬度值，将ncep四个网格点同一高程平面的大气参数进行空间上的线性插值；
(6)根据landsat8卫星过境的时间，将经过高程插值和空间插值后的大气参数在时间上进行线性插值，最终得到每个landsat8像元对应的大气参数。
本发明的有益效果是：针对山地的复杂地形，通过将坡度、坡向和天空可视因子等地形参数引入到热辐射传输方程，建立了山地地表温度反演方法，有效减小了利用常规的热辐射传输方程反演山地地表温度导致的误差，提高了山地地表温度反演精度。对比本发明与常规热辐射传输方程的反演结果，突出本发明在地形起伏较大的区域有明显的优势，能够弥补现有技术反演地表温度时造成的较大误差。图4表明，现有技术反演的地表温度在山谷处可造成1k的误差，这也进一步证实了本发明的有效效果。
附图说明
图1.山区地表温度反演流程图；
图2.山区地表温度反演结果图；
图3.山区地表高程图；
图4.山区地表温度反演结果与常规热辐射传输方程反演结果的温差图；
具体实施方式
以下结合具体实施例，对本发明进行详细说明。
步骤1：数据预处理
下载搭载在landsat8卫星上的tirs(thermalinfraredsensor)传感器的星上辐亮度产品，以及oli(operationallandimager)传感器的地表反射率产品。数据下载网址：https://earthexplorer.usgs.gov/。数据预处理包括:
(5)将星上辐亮度产品中的计数值转换为星上辐亮度：
l＝gain×dn offset(39)
式中，l为星上辐亮度；dn为计数值；gain和offset为增益和偏移。landsat8第10通道的增益和偏移分别为0.00033420和0.1。
(6)将地表反射率产品中的计数值转换为地表反射率：
ρ＝a×dn b(40)
式中，ρ为地表反射率，dn为计数值，a和b为转换系数。landsat8第4通道和第5通道的系数a和b分别为0.0000275和-0.2。
(7)计算归一化植被指数ndvi：
式中，ρred和ρnir分别为红光波段和近红外波段的地表反射率。
(8)计算植被覆盖度：
式中，pv为植被覆盖度，ndvi为当前像元的归一化植被指数，ndvimin和ndvimax分别为ndvi最小值和最大值，ndvimin取值为0.05，ndvimax取值为0.85。
步骤2：地表发射率计算
利用asterged产品计算landsat8第10通道的地表发射率。该产品是利用2000年到2008年所有的晴空aster数据通过温度与发射率分离算法反演的地表温度和发射率平均计算得到，包括aster通道10-14五个热红外通道发射率的平均值和标准偏差，以及其他数据集。数据下载网址：https://search.earthdata.nasa.gov/。landsat8地表发射率计算过程包括：
(1)从asterged产品中提取出ndvi、第13和14通道的dn值。根据公式(5)将dn值转换为asterged的地表发射率：
ε＝c×dn d(43)
式中，ε为asterged的地表发射率，dn为计数值，c和d为转换系数。aster第13和14通道的转换系数c和d分别为0.001和0。
(2)假定像元的地表发射率可以简单地表示为植被与裸土比辐射率的加权值：
εi＝εvipv εsi(1-pv)(44)
式中，εi为通道i的地表发射率，εvi为通道i的植被组分发射率，εsi为通道i的裸土组分发射率，pv为植被覆盖度。
(3)基于asterged的地表发射率和公式(6)，计算aster第13和14通道的裸土组分发射率：
式中，上标a代表aster传感器，i代表aster第13和14通道，
(4)从aster地物波谱库中挑选41条土壤反射率曲线，分别利用aster第13和14通道以及landsat8第10通道的光谱响应函数与反射率波谱曲线进行卷积计算，得到每个通道卷积后的土壤发射率。以aster第13和14通道卷积后的土壤发射率作为自变量，以landsat8第10通道卷积后的土壤发射率作为因变量，拟合得到三个通道土壤发射率之间的回归方程：
式中，
(5)landsat8地表发射率通过土壤组分发射率和植被组分发射率加权计得到：
式中，
步骤3：地形参数计算
下载地表高程dem数据。数据下载网址：https://search.earthdata.nasa.gov/。地形参数计算过程包括：
(1)计算坡度和坡向
坡度和坡向采用三阶反距离平方权差分法计算：
式中，s和a分别为坡度和坡向，slopens为以坡面目标点高程为中心的3×3个像元内的南北走向趋势，slopewe为东西走向的趋势，分别通过公式(12)和(13)计算：
slopewe＝(hj 1,i-1 2hj 1,i hj 1,i 1-hj-1,i-1-2hj-1,i-hj-1,i 1)/(8×cellsize)(50)
slopens＝(hj-1,i 1 2hj,i 1 hj 1,i 1-hj-1,i-1-2hj,i-1-hj 1,i-1)/(8×cellsize)(51)
式中，h为像元的地表高程，i和j分别为像元在地表高程图像中的行号和列号，cellsize为像元的空间分辨率。
(2)计算天空可视因子
天空可视因子表示山地地表一个点在半球空间内对天空的可视程度，即该点接收到的大气下行辐射与其在平坦地表未被遮挡时接收到的来自半球空间内的大气下行辐射的比值，取值范围为0到1，可以通过公式(14)计算：
式中，s为坡度，a为坡向，φ为方位角，hφ为水平张角，它表示当前目标像元所在水平面的法线与当前目标像元和周围邻近像元之间连线的夹角。由于我们难以将方位角360°方向的每个方向的天空可视因子都计算出来进行积分，所以我们将方位角离散为均匀分布的16个方向计算天空可视因子，这16个方向分别为0°,30°,45°,60°,90°,120°,135°,150°,180°,210°,225°,240°,270°,300°,315°,和330°。
步骤4：大气参数计算
下载ncep(nationalcentersforenvironmentalprediction)大气廓线数据。该数据提供每天utc时间00、06、12和18点四个时刻、空间分辨率为0.5°的大气温湿廓线。数据下载网址：https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs。大气参数计算过程包括：
(1)根据landsat8卫星过境的utc时间，选取最近两个时刻的大气廓线数据；
(2)提取ncep大气廓线31个气压层的高度、温度、相对湿度数据；
(3)将大气廓线输入大气辐射传输模型modtran，计算每条ncep大气廓线对应的大气参数，包括大气透过率、大气上行辐射和大气下行辐射；
(4)根据像元的地表高程，将距离该像元最近的ncep四个网格点的大气参数进行线性插值，插值到同一高程上；
(5)根据像元中心的经纬度值，将ncep四个网格点同一高程平面的大气参数进行空间上的线性插值；
(6)根据landsat8卫星过境的时间，将经过高程插值和空间插值后的大气参数在时间上进行线性插值，最终得到每个landsat8像元对应的大气参数。
步骤5：山地地表温度反演
在热红外波谱范围内，常规的热辐射传输方程表示为：
ltoa＝εb(ts)τ (1-ε)ldownτ lup(53)
式中，ltoa为大气层顶接收到的辐亮度，ε为地表发射率，b为普朗克函数，ts为地表温度，τ为大气透过率，ldown为大气下行辐射，lup为大气上行辐射。
山地地形的复杂几何结构使常规的热辐射传输模式发生变化。一方面，地形遮挡了部分来自半球方向的大气辐射，导致地表接收的大气下行辐射减少；另一方面，邻近地形的辐射增加了传感器接收的总辐亮度。山地热辐射传输方程表示为：
ltoa＝εb(ts)τ (1-ε)ldownτvd (1-ε)ladjτ lup(54)
式中，vd为天空可视因子，ladj为邻近地形辐射，表示为：
式中，下标a和b分别代表目标像元和邻近像元，lb为邻近像元的辐亮度，ta和tb分别为目标像元和邻近像元的法线与目标像元和邻近像元的中心点连线的夹角，dsa和dsb分别为目标像元和邻近像元的坡面面积，rab为目标像元和邻近像元之间的距离，n为邻近像元b的总个数。
基于山地热辐射传输方程反演山地地表温度表示为：
式中，tsm为反演的山地地表温度，b-1为普朗克函数的逆运算。
从公式(18)可以看出，在其他地表和大气参数都已知的情况下，想要反演得到山地地表温度，必须首先要知道邻近像元的辐亮度。想要知道邻近像元的辐亮度，就必须要知道邻近像元的地表温度。然而，邻近像元的地表温度事先是不知道的。为了解决这个问题，我们将常规热辐射传输方程反演得到的地表温度作为邻近像元地表温度的初始值，通过迭代计算的方法得到最终的山地地表温度。基于常规的热辐射传输方程反演地表温度表示为：
式中，tsf为反演的平坦地表温度。
具体的山地地表温度反演过程包括：
(1)使用步骤1中的landsat8第10通道的辐亮度、步骤2中计算的地表发射率，步骤4中计算的大气参数(大气透过率、大气上行辐射和大气下行辐射)，通过公式(19)逐像元反演得到平坦地表温度tsf，并将此时的地表温度作为邻近像元地表温度的初始值。
(2)将邻近像元地表温度的初始值输入公式(17)得到邻近像元辐亮度，并将步骤2中计算的地表发射率，步骤3中计算的坡度、坡向和天空可视因子，步骤4中计算的大气参数(大气透过率、大气上行辐射和大气下行辐射)，以及步骤1中的landsat8第10通道的辐亮度，同时输入公式(18)，逐像元反演得到山地地表温度tsm，并将此时的地表温度作为邻近像元的地表温度，用于迭代计算邻近像元辐亮度。
(3)重复过程(2)，直到前后两次邻近像元地表温度的差值小于0.01k或者迭代次数大于3时输出结果，得到最终的山地地表温度。
图2为采用本发明方法的山区地表温度反演结果图，该图表明，地表温度的空间分布与图3地形的空间分布有显著的相关性，且地表高程与地表温度的分布趋势相反，这表明本发明的结果符合温度直减率，即海拔每垂直升高1km，温度一般下降6.5k。同时有效说明反演山地地表温度时，地形几何参数的影响，突出了本发明的重要作用。
本发明的有效效果是：通过将坡度、坡向和天空可视因子等地形参数引入常规热辐射传输方程，考虑地形对热辐射能量组成部分的影响，建立山地地表热辐射传输方程，反演高精度的山地地表温度。对比本发明与常规热辐射传输方程的反演结果，突出本发明在地形起伏较大的区域有明显的优势，能够弥补现有技术反演地表温度时造成的较大误差。图4表明，现有技术反演的地表温度在山谷处可造成1k的误差，这也进一步证实了本发明的有效效果。
应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。
1.一种山地地表温度遥感反演方法，其特征在于，包括以下步骤：
步骤1：下载搭载在landsat8卫星上的tirs传感器的星上辐亮度产品，以及oli传感器的地表反射率产品，并进行数据预处理；
步骤2：地表发射率计算；
步骤3：地形参数计算；
下载地表高程dem数据，地形参数计算过程包括：
(1)计算坡度和坡向
(2)计算天空可视因子
步骤4：大气参数计算；大气参数包括大气透过率、大气上行辐射和大气下行辐射；
步骤5：山地地表温度反演；
在热红外波谱范围内，常规的热辐射传输方程表示为：
ltoa＝εb(ts)τ (1-ε)ldownτ lup(1)
式中，ltoa为大气层顶接收到的辐亮度，ε为地表发射率，b为普朗克函数，ts为地表温度，τ为大气透过率，ldown为大气下行辐射，lup为大气上行辐射；
山地热辐射传输方程表示为：
ltoa＝εb(ts)τ (1-ε)ldownτvd (1-ε)ladjτ lup(2)
式中，vd为天空可视因子，ladj为邻近地形辐射，表示为：
式中，下标a和b分别代表目标像元和邻近像元，lb为邻近像元的辐亮度，ta和tb分别为目标像元和邻近像元的法线与目标像元和邻近像元的中心点连线的夹角，dsa和dsb分别为目标像元和邻近像元的坡面面积，rab为目标像元和邻近像元之间的距离，n为邻近像元b的总个数；
基于山地热辐射传输方程反演山地地表温度表示为：
式中，tsm为反演的山地地表温度，b-1为普朗克函数的逆运算；
基于常规的热辐射传输方程反演地表温度表示为：
式中，tsf为反演的平坦地表温度。
2.根据权利要求1所述的山地地表温度遥感反演方法，其特征在于，所述步骤1中数据预处理包括:
(1)将星上辐亮度产品中的计数值转换为星上辐亮度：
l＝gain×dn offset(6)
式中，l为星上辐亮度；dn为计数值；gain和offset为增益和偏移；landsat8第10通道的增益和偏移分别为0.00033420和0.1；
(2)将地表反射率产品中的计数值转换为地表反射率：
ρ＝a×dn b(7)
式中，ρ为地表反射率，dn为计数值，a和b为转换系数；landsat8第4通道和第5通道的系数a和b分别为0.0000275和-0.2；
(3)计算归一化植被指数ndvi：
式中，ρred和ρnir分别为红光波段和近红外波段的地表反射率；
(4)计算植被覆盖度：
式中，pv为植被覆盖度，ndvi为当前像元的归一化植被指数，ndvimin和ndvimax分别为ndvi最小值和最大值，ndvimin取值为0.05，ndvimax取值为0.85。
3.根据权利要求1所述的山地地表温度遥感反演方法，其特征在于，所述步骤2中，利用asterged产品计算landsat8第10通道的地表发射率，该产品是利用2000年到2008年所有的晴空aster数据通过温度与发射率分离算法反演的地表温度和发射率平均计算得到，包括aster通道10-14五个热红外通道发射率的平均值和标准偏差，以及其他数据集；landsat8地表发射率计算过程包括：
(1)从asterged产品中提取出ndvi、第13和14通道的dn值；根据公式(5)将dn值转换为asterged的地表发射率：
ε＝c×dn d(10)
式中，ε为asterged的地表发射率，dn为计数值，c和d为转换系数；aster第13和14通道的转换系数c和d分别为0.001和0；
(2)假定像元的地表发射率可以简单地表示为植被与裸土比辐射率的加权值：
εi＝εvipv εsi(1-pv)(11)
式中，εi为通道i的地表发射率，εvi为通道i的植被组分发射率，εsi为通道i的裸土组分发射率，pv为植被覆盖度；
(3)基于asterged的地表发射率和公式(6)，计算aster第13和14通道的裸土组分发射率：
式中，上标a代表aster传感器，i代表aster第13和14通道，
(4)从aster地物波谱库中挑选41条土壤反射率曲线，分别利用aster第13和14通道以及landsat8第10通道的光谱响应函数与反射率波谱曲线进行卷积计算，得到每个通道卷积后的土壤发射率；以aster第13和14通道卷积后的土壤发射率作为自变量，以landsat8第10通道卷积后的土壤发射率作为因变量，拟合得到三个通道土壤发射率之间的回归方程：
式中，
(5)landsat8地表发射率通过土壤组分发射率和植被组分发射率加权计得到：
式中，
4.根据权利要求1所述的山地地表温度遥感反演方法，其特征在于，所述步骤4中，下载ncep大气廓线数据；大气参数计算过程包括：
(1)根据landsat8卫星过境的utc时间，选取最近两个时刻的大气廓线数据；
(2)提取ncep大气廓线31个气压层的高度、温度、相对湿度数据；
(3)将大气廓线输入大气辐射传输模型modtran，计算每条ncep大气廓线对应的大气参数，包括大气透过率、大气上行辐射和大气下行辐射；
(4)根据像元的地表高程，将距离该像元最近的ncep四个网格点的大气参数进行线性插值，插值到同一高程上；
(5)根据像元中心的经纬度值，将ncep四个网格点同一高程平面的大气参数进行空间上的线性插值；
(6)根据landsat8卫星过境的时间，将经过高程插值和空间插值后的大气参数在时间上进行线性插值，最终得到每个landsat8像元对应的大气参数。
5.根据权利要求1所述的山地地表温度遥感反演方法，其特征在于，所述计算坡度和坡向的方法为：坡度和坡向采用三阶反距离平方权差分法计算：
式中，s和a分别为坡度和坡向，slopens为以坡面目标点高程为中心的3×3个像元内的南北走向趋势，slopewe为东西走向的趋势，分别通过公式(12)和(13)计算：
slopewe＝(hj 1,i-1 2hj 1,i hj 1,i 1-hj-1,i-1-2hj-1,i-hj-1,i 1)/(8×cellsize)(17)
slopens＝(hj-1,i 1 2hj,i 1 hj 1,i 1-hj-1,i-1-2hj,i-1-hj 1,i-1)/(8×cellsize)(18)
式中，h为像元的地表高程，i和j分别为像元在地表高程图像中的行号和列号，cellsize为像元的空间分辨率。
6.根据权利要求5所述的山地地表温度遥感反演方法，其特征在于，所述计算天空可视因子方法为：
天空可视因子表示山地地表一个点在半球空间内对天空的可视程度，即该点接收到的大气下行辐射与其在平坦地表未被遮挡时接收到的来自半球空间内的大气下行辐射的比值，取值范围为0到1，通过公式(14)计算：
式中，s为坡度，a为坡向，φ为方位角，hφ为水平张角，它表示当前目标像元所在水平面的法线与当前目标像元和周围邻近像元之间连线的夹角。
7.根据权利要求7所述的山地地表温度遥感反演方法，其特征在于，所述方位角离散为均匀分布的16个方向计算天空可视因子，这16个方向分别为0°,30°,45°,60°,90°,120°,135°,150°,180°,210°,225°,240°,270°,300°,315°,和330°。
8.根据权利要求6所述的山地地表温度遥感反演方法，其特征在于，所述步骤5中，具体的山地地表温度反演过程：
(1)使用步骤1中的landsat8第10通道的辐亮度、步骤2中计算的地表发射率，步骤4中计算的大气参数，通过公式(19)逐像元反演得到平坦地表温度tsf，并将此时的地表温度作为邻近像元地表温度的初始值；
(2)将邻近像元地表温度的初始值输入公式(17)得到邻近像元辐亮度，并将步骤2中计算的地表发射率，步骤3中计算的坡度、坡向和天空可视因子，步骤4中计算的大气参数，以及步骤1中的landsat8第10通道的辐亮度，同时输入公式(18)，逐像元反演得到山地地表温度tsm，并将此时的地表温度作为邻近像元的地表温度，用于迭代计算邻近像元辐亮度；
(3)重复(2)，直到前后两次邻近像元地表温度的差值小于0.01k或者迭代次数大于3时输出结果，得到最终的山地地表温度。
技术总结