一种考虑城市三维几何结构影响的地表温度反演方法与流程

1.本发明涉及一种考虑城市三维几何结构影响的地表温度反演方法，属于定量遥感技术领域。


背景技术：

2.地表温度是地表与大气之间相互作用的关键因子，也是研究地表能量平衡的重要参数。了解地表温度的时空分布不仅在区域和全球尺度的地表能量过程中发挥着重要作用，而且对城市气候研究和城市热环境监测具有重要意义。卫星遥感技术为大范围获取高时空分辨率的地表温度提供了可能。目前，国内外学者已经发展了多种卫星数据地表温度遥感反演方法。然而，这些方法大多针对平坦地表，而对具有复杂几何结构的地表研究甚少。城市地表温度与人类活动息息相关，城市土地覆盖类型具有很强的异质性。不同类型的建筑物、人造地物和植被的复杂形态使得城市地表温度呈现出强烈的空间异质性。
3.受城市复杂三维几何结构的影响，城市地表温度的反演面临三个挑战：(1)像元内部的多次散射和反射对城市像元发射率的影响；(2)城市三维几何结构对大气向下辐射的影响；(3)邻近像元辐射对目标像元的影响。随着卫星传感器的空间分辨率越来越高，城市复杂三维几何结构对高空间分辨率热红外数据反演地表温度的影响越来越突出。
4.现有的地表温度反演方法大多针对平坦地表，主要适用于裸土、植被等自然地表，没有考虑到复杂下垫面以及地表几何结构对地表温度反演的影响。城市的复杂三维几何结构导致地表与大气之间的热辐射传输模式发生变化，对像元发射率、大气下行辐射以及邻近像元辐射等产生影响，使得基于常规热辐射传输方程的地表温度反演方法不能适用于具有复杂三维几何结构的城市地表，导致城市地表温度反演精度低。因此，现有技术存在缺陷，需要改进。


技术实现要素：

5.本发明所要解决的技术问题是针对现有技术的不足提供一种考虑城市三维几何结构影响的地表温度反演方法。
6.本发明的技术方案如下：
7.一种考虑城市三维几何结构影响的地表温度反演方法，包括以下步骤：城市地表温度反演过程包括：
8.步骤1：使用获得的landsat 8第10通道的辐亮度、计算得到的地表发射率以及计算得到的大气透过率、大气上行辐射和大气下行辐射，通过公式(10)逐像元反演得到平坦地表温度t
s_initial
，并将此时的地表温度作为邻近像元地表温度的初始值；
9.公式(10)：
[0010][0011]
式中，t
s_initial
为反演得到的平坦地表温度；
[0012]
步骤2：将邻近像元地表温度的初始值输入公式(8)得到邻近像元辐亮度，并将计算得到的天空可视因子，计算得到的城市地表的像元有效发射率，计算得到的大气透过率、大气上行辐射和大气下行辐射，以及获得的landsat 8第10通道的辐亮度，同时输入公式(9)，逐像元反演得到城市地表温度ts，并将此时的地表温度作为邻近像元的地表温度，用于迭代计算邻近像元辐亮度；
[0013]
公式(8)：
[0014][0015]
式中，下标a和b分别代表目标像元和邻近像元，lb为邻近像元的辐亮度，θ为目标像元和邻近像元的中心点连线的夹角，dsb为邻近像元的面积，r
ab
为目标像元和邻近像元之间的距离，n为邻近像元的总个数；
[0016]
公式(9)：
[0017][0018]
式中，ts为反演得到的考虑城市三维几何结构影响的地表温度，b-1
为普朗克函数的逆运算；
[0019]
步骤3：重复用于迭代计算邻近像元地表温度的步骤2，直到前后两次邻近像元地表温度的差值小于一定数值或者迭代次数大于一定数值时输出结果，得到最终的城市地表温度。
[0020]
所述的地表温度反演方法，所述landsat 8第10通道的辐亮度获取方法为：下载搭载在landsat 8卫星上的tirs传感器的星上辐亮度数据；数据下载网址：https://earthexplorer.usgs.gov/；利用公式(1)将星上辐亮度数据转换为星上辐亮度值：
[0021]
l＝gain
×
dn offset
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0022]
式中，l为星上辐亮度，dn为计数值，gain和offset分别为增益和偏移。对于landsat8第10通道，增益和偏移分别为0.00033420和0.1。
[0023]
所述的地表温度反演方法，所述天空可视因子的计算方法为：天空可视因子表示地表在半球空间内对天空的可视程度，即该点接收到的大气下行辐射与其在平坦未被遮挡时接收到的来自半球空间内的大气下行辐射的比值，取值范围为0到1，通过公式(2)计算得到：
[0024][0025]
式中，γi是给定搜索方向上的垂直仰角，n是搜索方向的数量。
[0026]
所述的地表温度反演方法，将整个半球空间的方位角离散为均匀分布的16个方向来计算天空可视因子，分别为0
°
,30
°
,45
°
,60
°
,90
°
,120
°
,135
°
,150
°
,180
°
,210
°
,225
°
,240
°
,270
°
,300
°
,315
°
和330
°
。
[0027]
所述的地表温度反演方法，所述地表发射率的计算方法为：
[0028]
(1)不考虑城市三维几何结构时，根据不同地物组分的材料发射率的加权平均值
计算得到平坦地表的像元发射率：
[0029][0030]
式中，ε
ef
为平坦地表的像元发射率，fk是像元内第k个地物组分的面积比例，εk是像元内第k个地物组分的材料发射率；
[0031]
(2)对于具有复杂三维几何结构的城市地表，像元有效发射率通过平坦地表的像元发射率ε
ef
与像元内部不同地物组分之间的多次散射和反射引起的发射率增量计算得到：
[0032][0033]
式中，ε
eu
为城市地表的像元有效发射率，δεm为像元内部不同地物组分之间的多次散射和反射引起的发射率增量；
[0034]
考虑像元内部不同地物组分之间的多次散射和反射，最终得到城市地表的像元有效发射率：
[0035][0036]
式中，υ
in
为landsat 8像元内部的天空可视因子，1-υ
in
为像元内被建筑物遮挡的可视天空的部分；当像元的天空可视因子接近0时，意味着天空对整个像元几乎是不可见的，辐射在像元内被无数次反射，像元的发射辐射趋于0，此时像元近似于一个黑体，像元的有效发射率为1；当像元的天空可视因子等于1时，意味着像元是一个没有任何建筑物遮挡的平坦地表，像元的有效发射率等于材料发射率。
[0037]
所述的地表温度反演方法，计算所述地表发射率之前，还包括城市地表覆盖类型图生成步骤：
[0038]
利用空间分辨率为3m的高空间分辨率卫星影像数据以及建筑物和水体矢量数据生成空间分辨率为3m的土地覆盖类型图；考虑到城市地表特征，将城市土地覆盖类型分为建筑物、道路、草地、树木和水体五类；高空间分辨率卫星影像和建筑物矢量数据从arctiler平台下载得到；高空间分辨率卫星图像来自quickbird商业卫星；水体矢量数据是从开放街道地图获得；将建筑物和水体矢量数据叠加在高空间分辨率卫星影像上，并使用基于高空间分辨率卫星影像的监督分类方法进一步区分道路、草地和树木，生成的城市土地覆盖类型图用于计算城市像元的有效发射率。
[0039]
所述的地表温度反演方法，所述大气透过率、大气上行辐射和大气下行辐射的计算方法为：将大气廓线数据输入到大气辐射传输软件modtran计算大气上行辐射、大气下行辐射和大气透过率等大气参数；
[0040]
大气参数逐像元计算过程包括：
[0041]
(1)根据landsat 8卫星的过境时间，选取最近时刻的大气廓线数据；
[0042]
(2)提取era5大气廓线37个气压层的高度、温度、相对湿度数据，并将大气廓线输入大气辐射传输模型modtran中，计算每条大气廓线对应的大气透过率、大气上行辐射和大气下行辐射三个参数；
[0043]
(3)根据像元的地表高程和中心经纬度值，将距离该像元最近的era5四个网格点的大气参数进行线性插值，插值到同一高程上，然后对同一高程平面的大气参数进行空间上的线性插值；
[0044]
(4)根据landsat 8卫星的过境时间，将经过高程和空间插值后的大气参数在时间上进行线性插值，最终得到每个landsat 8像元对应的大气参数。
[0045]
本发明的有效效果是：通过将天空可视因子等城市形态参数引入到常规热辐射传输方程，考虑城市三维几何结构对地表热辐射传输的影响，建立适用于城市地表的热辐射传输方程，反演高精度的城市地表温度。对比本发明与常规热辐射传输方程的地表温度反演结果，突出本发明在具有复杂三维几何结构的城市地表有明显的优势，能够弥补现有技术反演地表温度时造成的较大误差。图4表明，现有技术反演的地表温度在城市区域可造成1k的误差，这也进一步证实了本发明的有效效果。
附图说明
[0046]
图1为考虑城市三维几何结构影响的地表温度反演方法流程图；
[0047]
图2为研究区数字表面模型图；
[0048]
图3为考虑城市三维几何结构影响的地表温度反演结果图，(a)2019年8月17日，(b)2019年12月7日；
[0049]
图4为考虑城市三维几何结构影响的地表温度反演结果与基于常规热辐射传输方程的地表温度反演结果之间的温差图，(a)2019年8月17日，(b)2019年12月7日；
具体实施方式
[0050]
以下结合具体实施例，对本发明进行详细说明。
[0051]
步骤1：星上辐亮度产品获取和数据预处理
[0052]
下载搭载在landsat 8卫星上的tirs(thermal infrared sensor)传感器的星上辐亮度数据。数据下载网址：https://earthexplorer.usgs.gov/。利用公式(1)将星上辐亮度数据转换为星上辐亮度值：
[0053]
l＝gain
×
dn offset
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0054]
式中，l为星上辐亮度，dn为计数值，gain和offset分别为增益和偏移。对于landsat 8第10通道，增益和偏移分别为0.00033420和0.1。
[0055]
步骤2：天空可视因子计算
[0056]
利用数字表面模型(digital surface model,dsm)数据计算天空可视因子(sky view factor,svf)。从自然资源卫星遥感云服务平台(http://www.sasclouds.com/en)获得空间分辨率为3m的dsm数据。该数据由国产资源三号卫星立体影像生成，可以提供建筑物和植物的高度信息。
[0057]
天空可视因子表示地表在半球空间内对天空的可视程度，即该点接收到的大气下行辐射与其在平坦未被遮挡时接收到的来自半球空间内的大气下行辐射的比值，取值范围为0到1，通过公式(2)计算得到：
[0058]
[0059]
式中，γi是给定搜索方向上的垂直仰角，n是搜索方向的数量。
[0060]
将整个半球空间的方位角离散为均匀分布的16个方向来计算天空可视因子，分别为0
°
,30
°
,45
°
,60
°
,90
°
,120
°
,135
°
,150
°
,180
°
,210
°
,225
°
,240
°
,270
°
,300
°
,315
°
和330
°
。
[0061]
步骤3：城市地表覆盖类型图生成
[0062]
利用空间分辨率为3m的高空间分辨率卫星影像数据以及建筑物和水体矢量数据生成空间分辨率为3m的土地覆盖类型图。考虑到城市地表特征，将城市土地覆盖类型分为建筑物、道路、草地、树木和水体五类。高空间分辨率卫星影像和建筑物矢量数据从arctiler平台下载得到(http://www.arctiler.com/)。高空间分辨率卫星图像来自quickbird商业卫星。水体矢量数据是从开放街道地图(https://www.openstreetmap.org/)获得。将建筑物和水体矢量数据叠加在高空间分辨率卫星影像上，并使用基于高空间分辨率卫星影像的监督分类方法进一步区分道路、草地和树木，生成的城市土地覆盖类型图用于计算城市像元的有效发射率。
[0063]
步骤4：城市像元有效发射率计算
[0064]
利用aster波谱库和城市不透水材料波谱库计算城市典型地物的材料发射率。除了自然地表(草地、树木和水体)，城市地表主要被不透水面覆盖，例如由沥青或混凝土构成的道路，由混凝土和玻璃构成的建筑物等。为了获得城市典型地物的材料发射率，将建筑物的材料发射率设置为砖石、混凝土和玻璃的发射率平均值，将道路的材料发射率设置为水泥、混凝土和砖石的发射率平均值。砖石、水泥、混凝土和玻璃的材料发射率通过landsat 8的传感器响应函数与城市不透水材料波谱库中的砖石、水泥、混凝土和玻璃的发射光谱曲线卷积计算得到。自然地表中草地、树木和水体的材料发射率通过landsat 8的传感器响应函数与aster波谱库中的草地、树木和水体的发射光谱曲线卷积计算得到。对于平坦的地表，像元发射率由像元内部不同地物组分的材料发射率的加权平均值计算得到。对于具有复杂三维几何结构的城市地表，则需要考虑像元内部不同地物组分之间的多次反射和散射计算城市地表的像元有效发射率。城市地表的像元有效发射率计算过程包括：
[0065]
(1)不考虑城市三维几何结构时，根据不同地物组分的材料发射率的加权平均值计算得到平坦地表的像元发射率：
[0066][0067]
式中，ε
ef
为平坦地表的像元发射率，fk是像元内第k个地物组分的面积比例，εk是像元内第k个地物组分的材料发射率。
[0068]
(2)对于具有复杂三维几何结构的城市地表，像元有效发射率通过平坦地表的像元发射率(即公式(3)中的ε
ef
)与像元内部不同地物组分之间的多次散射和反射(即腔体效应)引起的发射率增量计算得到：
[0069][0070]
式中，ε
eu
为城市地表的像元有效发射率，δεm为像元内部不同地物组分之间的多次散射和反射引起的发射率增量。
[0071]
考虑像元内部不同地物组分之间的多次散射和反射，最终得到城市地表的像元有
效发射率：
[0072][0073]
式中，υ
in
为landsat 8像元内部的天空可视因子，1-υ
in
为像元内被建筑物遮挡的可视天空的部分。当像元的天空可视因子接近0时，意味着天空对整个像元几乎是不可见的，辐射在像元内被无数次反射，像元的发射辐射趋于0，此时像元近似于一个黑体，像元的有效发射率为1；当像元的天空可视因子等于1时，意味着像元是一个没有任何建筑物遮挡的平坦地表，像元的有效发射率等于材料发射率。
[0074]
步骤5：大气参数逐像元计算
[0075]
era5再分析廓线数据是由欧洲中期天气预报中心(european centre medium-range weather forecasts,ecmwf)发布的最新一代的再分析数据。era5在其前一代era-interim的基础上实现了时间和空间分辨率的提升。该数据的空间分辨率为0.25
°×
0.25
°
，时间分辨率为1小时。era5提供1hpa至1000hpa的37个压强层的大气温度和湿度廓线数据。将大气廓线数据输入到大气辐射传输软件modtran计算大气上行辐射、大气下行辐射和大气透过率等大气参数。
[0076]
大气参数逐像元计算过程包括：
[0077]
(1)根据landsat 8卫星的过境时间，选取最近时刻的大气廓线数据；
[0078]
(2)提取era5大气廓线37个气压层的高度、温度、相对湿度数据，并将大气廓线输入大气辐射传输模型modtran中，计算每条大气廓线对应的大气透过率、大气上行辐射和大气下行辐射三个参数；
[0079]
(3)根据像元的地表高程和中心经纬度值，将距离该像元最近的era5四个网格点的大气参数进行线性插值，插值到同一高程上，然后对同一高程平面的大气参数进行空间上的线性插值；
[0080]
(4)根据landsat 8卫星的过境时间，将经过高程和空间插值后的大气参数在时间上进行线性插值，最终得到每个landsat 8像元对应的大气参数。
[0081]
步骤6：城市地表温度反演
[0082]
在热红外波谱范围内，常规的热辐射传输方程表示为：
[0083]
l
toa
＝[ε
ef
b(ts) (1-ε
ef
)ld]τ luꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0084]
式中，l
toa
为大气层顶接收到的辐亮度，ε
ef
为平坦地表的像元发射率，b为普朗克函数，ts为地表温度，τ为大气透过率，ld为大气下行辐射，lu为大气上行辐射。
[0085]
城市地表的复杂三维几何结构使常规的热辐射传输模式发生变化。一方面，邻近像元遮挡了部分来自半球方向的大气辐射，导致地表接收的大气下行辐射减少；另一方面，来自邻近像元的辐射增加了传感器接收的总辐射。因此，考虑城市三维几何结构影响的热辐射传输方程表示为：
[0086]
l
toa
＝[ε
eu
b(ts) (1-ε
eu
)(υ
out
ld l
adj
)]τ luꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0087]
式中，l
toa
为大气层顶接收到的辐亮度，ε
eu
为城市地表的像元有效发射率，υ
out
为landsat8像元之间的天空可视因子，l
adj
为邻近像元辐射，表示为：
[0088][0089]
式中，下标a和b分别代表目标像元和邻近像元，lb为邻近像元的辐亮度，θ为目标像元和邻近像元的中心点连线的夹角，dsb为邻近像元的面积，r
ab
为目标像元和邻近像元之间的距离，n为邻近像元的总个数。
[0090]
基于考虑城市三维几何结构影响的热辐射传输方程反演的城市地表温度表示为：
[0091][0092]
式中，ts为反演得到的考虑城市三维几何结构影响的地表温度，b-1
为普朗克函数的逆运算。
[0093]
从公式(9)可以看出，在其他地表和大气参数都已知的情况下，首先要获取邻近像元的辐亮度。计算邻近像元辐亮度的前提是获取邻近像元的地表温度。然而，邻近像元的地表温度事先是未知的。为了解决这一问题，将常规热辐射传输方程反演得到的地表温度作为邻近像元地表温度的初始值，通过迭代计算的方法得到最终的城市地表温度。基于常规的热辐射传输方程反演的地表温度表示为：
[0094][0095]
式中，t
s_initial
为反演得到的平坦地表温度。
[0096]
城市地表温度反演过程包括：
[0097]
(1)使用步骤1中的landsat 8第10通道的辐亮度、步骤4中计算的地表发射率以及步骤5中计算的大气透过率、大气上行辐射和大气下行辐射，通过公式(10)逐像元反演得到平坦地表温度t
s_initial
，并将此时的地表温度作为邻近像元地表温度的初始值。
[0098]
(2)将邻近像元地表温度的初始值输入公式(8)得到邻近像元辐亮度，并将步骤2中计算得到的天空可视因子，步骤4中计算的城市地表的像元有效发射率，步骤5中计算的大气透过率、大气上行辐射和大气下行辐射，以及步骤1中的landsat 8第10通道的辐亮度，同时输入公式(9)，逐像元反演得到城市地表温度ts，并将此时的地表温度作为邻近像元的地表温度，用于迭代计算邻近像元辐亮度。
[0099]
(3)重复用于迭代计算邻近像元地表温度的过程(2)，直到前后两次邻近像元地表温度的差值小于0.1k或者迭代次数大于3时输出结果，得到最终的城市地表温度。
[0100]
图3为采用本发明方法反演的城市地表温度结果图。该图表明，地表温度的空间分布与建筑物高度和地表覆盖类型有关，其中，建筑物高度较低的区域的地表温度较高，水体及植被覆盖区域的地表温度较低。图4表明，基于常规热辐射传输方程的地表温度反演方法会对地表温度的反演结果造成最高1k的误差，同时有效说明在精确反演城市地表温度时，考虑城市三维几何结构影响的必要性，突出了本发明的重要作用。
[0101]
应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。