CN112487346A - 一种山地地表温度遥感反演方法 - Google Patents
一种山地地表温度遥感反演方法 Download PDFInfo
- Publication number
- CN112487346A CN112487346A CN202011153753.8A CN202011153753A CN112487346A CN 112487346 A CN112487346 A CN 112487346A CN 202011153753 A CN202011153753 A CN 202011153753A CN 112487346 A CN112487346 A CN 112487346A
- Authority
- CN
- China
- Prior art keywords
- emissivity
- surface temperature
- atmospheric
- aster
- channel
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 44
- 230000005855 radiation Effects 0.000 claims abstract description 63
- 230000005540 biological transmission Effects 0.000 claims abstract description 21
- 238000002310 reflectometry Methods 0.000 claims abstract description 14
- 238000002834 transmittance Methods 0.000 claims abstract description 8
- 238000007781 pre-processing Methods 0.000 claims abstract description 6
- 238000000179 transient infrared spectroscopy Methods 0.000 claims abstract description 3
- 241000132092 Aster Species 0.000 claims description 69
- 239000002689 soil Substances 0.000 claims description 34
- 230000003595 spectral effect Effects 0.000 claims description 24
- 238000004364 calculation method Methods 0.000 claims description 21
- 230000008569 process Effects 0.000 claims description 13
- 238000005316 response function Methods 0.000 claims description 9
- 230000006870 function Effects 0.000 claims description 7
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 150000001875 compounds Chemical class 0.000 claims description 6
- 238000001228 spectrum Methods 0.000 claims description 6
- 230000001419 dependent effect Effects 0.000 claims description 3
- 230000035699 permeability Effects 0.000 claims description 3
- 238000000926 separation method Methods 0.000 claims description 3
- 239000006185 dispersion Substances 0.000 claims 1
- 230000008859 change Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
- G01J5/007—Radiation pyrometry, e.g. infrared or optical thermometry for earth observation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
- G01J5/80—Calibration
- G01J5/802—Calibration by correcting for emissivity
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
- G01J5/80—Calibration
- G01J5/804—Calibration using atmospheric correction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01K—MEASURING TEMPERATURE; MEASURING QUANTITY OF HEAT; THERMALLY-SENSITIVE ELEMENTS NOT OTHERWISE PROVIDED FOR
- G01K11/00—Measuring temperature based upon physical or chemical changes not covered by groups G01K3/00, G01K5/00, G01K7/00 or G01K9/00
- G01K11/12—Measuring temperature based upon physical or chemical changes not covered by groups G01K3/00, G01K5/00, G01K7/00 or G01K9/00 using changes in colour, translucency or reflectance
- G01K11/125—Measuring temperature based upon physical or chemical changes not covered by groups G01K3/00, G01K5/00, G01K7/00 or G01K9/00 using changes in colour, translucency or reflectance using changes in reflectance
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Operations Research (AREA)
- Life Sciences & Earth Sciences (AREA)
- Geology (AREA)
- Probability & Statistics with Applications (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Computing Systems (AREA)
- Radiation Pyrometers (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种山地地表温度遥感反演方法,包括以下步骤:步骤1:下载搭载在Landsat 8卫星上的TIRS传感器的星上辐亮度产品,以及OLI传感器的地表反射率产品,并进行数据预处理;步骤2:地表发射率计算;步骤3:地形参数计算;步骤4:大气参数计算;大气参数包括大气透过率、大气上行辐射和大气下行辐射;步骤5:山地地表温度反演;针对山地的复杂地形,通过将坡度、坡向和天空可视因子等地形参数引入到热辐射传输方程,建立了山地地表温度反演方法,有效减小了利用常规的热辐射传输方程反演山地地表温度导致的误差,提高了山地地表温度反演精度。
Description
技术领域
本发明涉及遥感技术领域,尤其涉及的是一种山地地表温度遥感反演方法。
背景技术
地表温度是描述地表能量分布情况以及地表与大气相互作用的关键因素,也是研究地表物质能量动态平衡的重要参数,被广泛应用于灾害评估、蒸散发估算、气候研究和土壤水分监测等方面。目前的地表温度反演方法大多针对平坦地面,开展复杂地形区的研究甚少。我国五大典型地形包括平原、山地、高原、盆地和丘陵,其中平原面积仅占12%。在复杂地形区,由于地形的起伏造成地表结构发生变化,这种变化引起地面与传感器的几何关系发生改变,进而导致传感器接收到的地面辐射能量的组成部分发生变化。现如今,越来越多高空间分辨率的卫星数据盛行,地形影响在高空间分辨率数据中尤为突出,这使得我们不得不考虑地形起伏导致的辐射能量组成部分的变化。
现有地表温度反演方法主要适用于平坦地表,针对复杂地形区的研究甚少。然而,山地地形的复杂几何结构使得地表与大气之间的热辐射传输方式发生变化,这使得基于常规热辐射传输方程的地表温度反演方法不能适用于山地地表温度反演。
因此,现有技术存在缺陷,需要改进。
发明内容
本发明所要解决的技术问题是针对现有技术的不足提供一种山地地表温度遥感反演方法。
本发明的技术方案如下:
一种山地地表温度遥感反演方法,包括以下步骤:
步骤1:下载搭载在Landsat 8卫星上的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中的Landsat 8第10通道的辐亮度、步骤2中计算的地表发射率,步骤4中计算的大气参数,通过公式(19)逐像元反演得到平坦地表温度Tsf,并将此时的地表温度作为邻近像元地表温度的初始值;
(2)将邻近像元地表温度的初始值输入公式(17)得到邻近像元辐亮度,并将步骤2中计算的地表发射率,步骤3中计算的坡度、坡向和天空可视因子,步骤4中计算的大气参数,以及步骤1中的Landsat 8第10通道的辐亮度,同时输入公式(18),逐像元反演得到山地地表温度Tsm,并将此时的地表温度作为邻近像元的地表温度,用于迭代计算邻近像元辐亮度;
(3)重复(2),直到前后两次邻近像元地表温度的差值小于0.01K或者迭代次数大于3时输出结果,得到最终的山地地表温度。
所述的山地地表温度遥感反演方法,所述步骤1中数据预处理包括:
(1)将星上辐亮度产品中的计数值转换为星上辐亮度:
L=gain×DN+offset (30)
式中,L为星上辐亮度;DN为计数值;gain和offset为增益和偏移;Landsat 8第10通道的增益和偏移分别为0.00033420和0.1;
(2)将地表反射率产品中的计数值转换为地表反射率:
ρ=a×DN+b (31)
式中,ρ为地表反射率,DN为计数值,a和b为转换系数;Landsat 8第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中,利用ASTER GED产品计算Landsat 8第10通道的地表发射率,该产品是利用2000年到2008年所有的晴空ASTER数据通过温度与发射率分离算法反演的地表温度和发射率平均计算得到,包括ASTER通道10-14五个热红外通道发射率的平均值和标准偏差,以及其他数据集;Landsat 8地表发射率计算过程包括:
(1)从ASTER GED产品中提取出NDVI、第13和14通道的DN值;根据公式(5)将DN值转换为ASTER GED的地表发射率:
ε=c×DN+d (34)
式中,ε为ASTER GED的地表发射率,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)基于ASTER GED的地表发射率和公式(6),计算ASTER第13和14通道的裸土组分发射率:
式中,上标A代表ASTER传感器,i代表ASTER第13和14通道,为通道i的ASTER GED的地表发射率,为通道i的ASTER的植被组分发射率,可通过ASTER的光谱响应函数与ASTER地物波谱库的植被波谱曲线进行卷积计算得到,为通道i的ASTER的裸土组分发射率,为ASTER的植被覆盖度,可基于ASTER GED的NDVI通过公式(4)计算得到;
(4)从ASTER地物波谱库中挑选41条土壤反射率曲线,分别利用ASTER第13和14通道以及Landsat 8第10通道的光谱响应函数与反射率波谱曲线进行卷积计算,得到每个通道卷积后的土壤发射率;以ASTER第13和14通道卷积后的土壤发射率作为自变量,以Landsat 8第10通道卷积后的土壤发射率作为因变量,拟合得到三个通道土壤发射率之间的回归方程:
(5)Landsat 8地表发射率通过土壤组分发射率和植被组分发射率加权计得到:
式中,为Landsat 8第10通道的地表发射率,为Landsat 8第10通道的土壤组分发射率,为Landsat 8第10通道的植被组分发射率,可通过Landsat 8第10通道的光谱响应函数与ASTER地物波谱库的植被波谱曲线进行卷积计算得到,为Landsat 8的植被覆盖度,可基于Landsat 8的地表发射率通过公式(3)计算NDVI,并进一步通过公式(4)得到。
所述的山地地表温度遥感反演方法,所述步骤4中,下载NCEP大气廓线数据;大气参数计算过程包括:
(1)根据Landsat 8卫星过境的UTC时间,选取最近两个时刻的大气廓线数据;
(2)提取NCEP大气廓线31个气压层的高度、温度、相对湿度数据;
(3)将大气廓线输入大气辐射传输模型MODTRAN,计算每条NCEP大气廓线对应的大气参数,包括大气透过率、大气上行辐射和大气下行辐射;
(4)根据像元的地表高程,将距离该像元最近的NCEP四个网格点的大气参数进行线性插值,插值到同一高程上;
(5)根据像元中心的经纬度值,将NCEP四个网格点同一高程平面的大气参数进行空间上的线性插值;
(6)根据Landsat 8卫星过境的时间,将经过高程插值和空间插值后的大气参数在时间上进行线性插值,最终得到每个Landsat 8像元对应的大气参数。
本发明的有益效果是:针对山地的复杂地形,通过将坡度、坡向和天空可视因子等地形参数引入到热辐射传输方程,建立了山地地表温度反演方法,有效减小了利用常规的热辐射传输方程反演山地地表温度导致的误差,提高了山地地表温度反演精度。对比本发明与常规热辐射传输方程的反演结果,突出本发明在地形起伏较大的区域有明显的优势,能够弥补现有技术反演地表温度时造成的较大误差。图4表明,现有技术反演的地表温度在山谷处可造成1K的误差,这也进一步证实了本发明的有效效果。
附图说明
图1.山区地表温度反演流程图;
图2.山区地表温度反演结果图;
图3.山区地表高程图;
图4.山区地表温度反演结果与常规热辐射传输方程反演结果的温差图;
具体实施方式
以下结合具体实施例,对本发明进行详细说明。
步骤1:数据预处理
下载搭载在Landsat 8卫星上的TIRS(Thermal Infrared Sensor)传感器的星上辐亮度产品,以及OLI(Operational Land Imager)传感器的地表反射率产品。数据下载网址:https://earthexplorer.usgs.gov/。数据预处理包括:
(5)将星上辐亮度产品中的计数值转换为星上辐亮度:
L=gain×DN+offset (39)
式中,L为星上辐亮度;DN为计数值;gain和offset为增益和偏移。Landsat 8第10通道的增益和偏移分别为0.00033420和0.1。
(6)将地表反射率产品中的计数值转换为地表反射率:
ρ=a×DN+b (40)
式中,ρ为地表反射率,DN为计数值,a和b为转换系数。Landsat 8第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:地表发射率计算
利用ASTER GED产品计算Landsat 8第10通道的地表发射率。该产品是利用2000年到2008年所有的晴空ASTER数据通过温度与发射率分离算法反演的地表温度和发射率平均计算得到,包括ASTER通道10-14五个热红外通道发射率的平均值和标准偏差,以及其他数据集。数据下载网址:https://search.earthdata.nasa.gov/。Landsat 8地表发射率计算过程包括:
(1)从ASTER GED产品中提取出NDVI、第13和14通道的DN值。根据公式(5)将DN值转换为ASTER GED的地表发射率:
ε=c×DN+d (43)
式中,ε为ASTER GED的地表发射率,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)基于ASTER GED的地表发射率和公式(6),计算ASTER第13和14通道的裸土组分发射率:
式中,上标A代表ASTER传感器,i代表ASTER第13和14通道,为通道i的ASTER GED的地表发射率,为通道i的ASTER的植被组分发射率,可通过ASTER的光谱响应函数与ASTER地物波谱库的植被波谱曲线进行卷积计算得到,为通道i的ASTER的裸土组分发射率,为ASTER的植被覆盖度,可基于ASTER GED的NDVI通过公式(4)计算得到。
(4)从ASTER地物波谱库中挑选41条土壤反射率曲线,分别利用ASTER第13和14通道以及Landsat 8第10通道的光谱响应函数与反射率波谱曲线进行卷积计算,得到每个通道卷积后的土壤发射率。以ASTER第13和14通道卷积后的土壤发射率作为自变量,以Landsat 8第10通道卷积后的土壤发射率作为因变量,拟合得到三个通道土壤发射率之间的回归方程:
(5)Landsat 8地表发射率通过土壤组分发射率和植被组分发射率加权计得到:
式中,为Landsat 8第10通道的地表发射率,为Landsat 8第10通道的土壤组分发射率,为Landsat 8第10通道的植被组分发射率,可通过Landsat 8第10通道的光谱响应函数与ASTER地物波谱库的植被波谱曲线进行卷积计算得到,为Landsat 8的植被覆盖度,可基于Landsat 8的地表发射率通过公式(3)计算NDVI,并进一步通过公式(4)得到。
步骤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(National Centers for Environmental Prediction)大气廓线数据。该数据提供每天UTC时间00、06、12和18点四个时刻、空间分辨率为0.5°的大气温湿廓线。数据下载网址:https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forc ast-system-gfs。大气参数计算过程包括:
(1)根据Landsat 8卫星过境的UTC时间,选取最近两个时刻的大气廓线数据;
(2)提取NCEP大气廓线31个气压层的高度、温度、相对湿度数据;
(3)将大气廓线输入大气辐射传输模型MODTRAN,计算每条NCEP大气廓线对应的大气参数,包括大气透过率、大气上行辐射和大气下行辐射;
(4)根据像元的地表高程,将距离该像元最近的NCEP四个网格点的大气参数进行线性插值,插值到同一高程上;
(5)根据像元中心的经纬度值,将NCEP四个网格点同一高程平面的大气参数进行空间上的线性插值;
(6)根据Landsat 8卫星过境的时间,将经过高程插值和空间插值后的大气参数在时间上进行线性插值,最终得到每个Landsat 8像元对应的大气参数。
步骤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中的Landsat 8第10通道的辐亮度、步骤2中计算的地表发射率,步骤4中计算的大气参数(大气透过率、大气上行辐射和大气下行辐射),通过公式(19)逐像元反演得到平坦地表温度Tsf,并将此时的地表温度作为邻近像元地表温度的初始值。
(2)将邻近像元地表温度的初始值输入公式(17)得到邻近像元辐亮度,并将步骤2中计算的地表发射率,步骤3中计算的坡度、坡向和天空可视因子,步骤4中计算的大气参数(大气透过率、大气上行辐射和大气下行辐射),以及步骤1中的Landsat 8第10通道的辐亮度,同时输入公式(18),逐像元反演得到山地地表温度Tsm,并将此时的地表温度作为邻近像元的地表温度,用于迭代计算邻近像元辐亮度。
(3)重复过程(2),直到前后两次邻近像元地表温度的差值小于0.01K或者迭代次数大于3时输出结果,得到最终的山地地表温度。
图2为采用本发明方法的山区地表温度反演结果图,该图表明,地表温度的空间分布与图3地形的空间分布有显著的相关性,且地表高程与地表温度的分布趋势相反,这表明本发明的结果符合温度直减率,即海拔每垂直升高1km,温度一般下降6.5K。同时有效说明反演山地地表温度时,地形几何参数的影响,突出了本发明的重要作用。
本发明的有效效果是:通过将坡度、坡向和天空可视因子等地形参数引入常规热辐射传输方程,考虑地形对热辐射能量组成部分的影响,建立山地地表热辐射传输方程,反演高精度的山地地表温度。对比本发明与常规热辐射传输方程的反演结果,突出本发明在地形起伏较大的区域有明显的优势,能够弥补现有技术反演地表温度时造成的较大误差。图4表明,现有技术反演的地表温度在山谷处可造成1K的误差,这也进一步证实了本发明的有效效果。
应当理解的是,对本领域普通技术人员来说,可以根据上述说明加以改进或变换,而所有这些改进和变换都应属于本发明所附权利要求的保护范围。
Claims (8)
1.一种山地地表温度遥感反演方法,其特征在于,包括以下步骤:
步骤1:下载搭载在Landsat 8卫星上的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为增益和偏移;Landsat 8第10通道的增益和偏移分别为0.00033420和0.1;
(2)将地表反射率产品中的计数值转换为地表反射率:
ρ=a×DN+b (7)
式中,ρ为地表反射率,DN为计数值,a和b为转换系数;Landsat 8第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中,利用ASTER GED产品计算Landsat 8第10通道的地表发射率,该产品是利用2000年到2008年所有的晴空ASTER数据通过温度与发射率分离算法反演的地表温度和发射率平均计算得到,包括ASTER通道10-14五个热红外通道发射率的平均值和标准偏差,以及其他数据集;Landsat8地表发射率计算过程包括:
(1)从ASTER GED产品中提取出NDVI、第13和14通道的DN值;根据公式(5)将DN值转换为ASTER GED的地表发射率:
ε=c×DN+d (10)
式中,ε为ASTER GED的地表发射率,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)基于ASTER GED的地表发射率和公式(6),计算ASTER第13和14通道的裸土组分发射率:
式中,上标A代表ASTER传感器,i代表ASTER第13和14通道,为通道i的ASTER GED的地表发射率,为通道i的ASTER的植被组分发射率,可通过ASTER的光谱响应函数与ASTER地物波谱库的植被波谱曲线进行卷积计算得到,为通道i的ASTER的裸土组分发射率,为ASTER的植被覆盖度;
(4)从ASTER地物波谱库中挑选41条土壤反射率曲线,分别利用ASTER第13和14通道以及Landsat 8第10通道的光谱响应函数与反射率波谱曲线进行卷积计算,得到每个通道卷积后的土壤发射率;以ASTER第13和14通道卷积后的土壤发射率作为自变量,以Landsat 8第10通道卷积后的土壤发射率作为因变量,拟合得到三个通道土壤发射率之间的回归方程:
(5)Landsat 8地表发射率通过土壤组分发射率和植被组分发射率加权计得到:
4.根据权利要求1所述的山地地表温度遥感反演方法,其特征在于,所述步骤4中,下载NCEP大气廓线数据;大气参数计算过程包括:
(1)根据Landsat 8卫星过境的UTC时间,选取最近两个时刻的大气廓线数据;
(2)提取NCEP大气廓线31个气压层的高度、温度、相对湿度数据;
(3)将大气廓线输入大气辐射传输模型MODTRAN,计算每条NCEP大气廓线对应的大气参数,包括大气透过率、大气上行辐射和大气下行辐射;
(4)根据像元的地表高程,将距离该像元最近的NCEP四个网格点的大气参数进行线性插值,插值到同一高程上;
(5)根据像元中心的经纬度值,将NCEP四个网格点同一高程平面的大气参数进行空间上的线性插值;
(6)根据Landsat 8卫星过境的时间,将经过高程插值和空间插值后的大气参数在时间上进行线性插值,最终得到每个Landsat 8像元对应的大气参数。
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为像元的空间分辨率。
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中的Landsat 8第10通道的辐亮度、步骤2中计算的地表发射率,步骤4中计算的大气参数,通过公式(19)逐像元反演得到平坦地表温度Tsf,并将此时的地表温度作为邻近像元地表温度的初始值;
(2)将邻近像元地表温度的初始值输入公式(17)得到邻近像元辐亮度,并将步骤2中计算的地表发射率,步骤3中计算的坡度、坡向和天空可视因子,步骤4中计算的大气参数,以及步骤1中的Landsat 8第10通道的辐亮度,同时输入公式(18),逐像元反演得到山地地表温度Tsm,并将此时的地表温度作为邻近像元的地表温度,用于迭代计算邻近像元辐亮度;
(3)重复(2),直到前后两次邻近像元地表温度的差值小于0.01K或者迭代次数大于3时输出结果,得到最终的山地地表温度。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011153753.8A CN112487346B (zh) | 2020-10-26 | 2020-10-26 | 一种山地地表温度遥感反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011153753.8A CN112487346B (zh) | 2020-10-26 | 2020-10-26 | 一种山地地表温度遥感反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112487346A true CN112487346A (zh) | 2021-03-12 |
CN112487346B CN112487346B (zh) | 2021-07-23 |
Family
ID=74927572
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011153753.8A Active CN112487346B (zh) | 2020-10-26 | 2020-10-26 | 一种山地地表温度遥感反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112487346B (zh) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113269464A (zh) * | 2021-06-10 | 2021-08-17 | 中国科学院地理科学与资源研究所 | 一种生态恢复评估方法和生态恢复评估装置 |
CN113449438A (zh) * | 2021-07-30 | 2021-09-28 | 北京环境特性研究所 | 一种可见光至热红外波段对地遥感图像仿真方法 |
CN113743000A (zh) * | 2021-08-13 | 2021-12-03 | 电子科技大学 | 一种生成高时间分辨率全天候地表温度的方法 |
CN113887780A (zh) * | 2021-08-26 | 2022-01-04 | 国家卫星气象中心(国家空间天气监测预警中心) | 一种卫星遥感地表温度的估计方法、装置及设备 |
CN113945282A (zh) * | 2021-09-03 | 2022-01-18 | 中国空间技术研究院 | 红外遥感卫星温度反演精度指标分配系统及方法 |
CN114152350A (zh) * | 2021-12-09 | 2022-03-08 | 中国农业科学院农业资源与农业区划研究所 | 一种考虑城市三维几何结构影响的地表温度反演方法 |
CN115374390A (zh) * | 2022-10-26 | 2022-11-22 | 中国科学院地理科学与资源研究所 | 一种高光谱热红外地表温度和发射率波谱协同反演方法 |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR101207925B1 (ko) * | 2010-02-22 | 2012-12-04 | 공주대학교 산학협력단 | 지표면 온도 산출 방법 |
WO2013077557A1 (ko) * | 2011-11-21 | 2013-05-30 | 한국항공우주연구원 | 성층권 비행선과 반사판을 이용한 지표면의 온도 조절 방법 |
CN103810387A (zh) * | 2014-02-13 | 2014-05-21 | 中国科学院地理科学与资源研究所 | 基于modis数据的地表蒸散发全遥感反演方法及系统 |
CN104657935A (zh) * | 2015-03-09 | 2015-05-27 | 广西壮族自治区气象减灾研究所 | 一种近地面气温反演方法 |
CN109509319A (zh) * | 2018-12-29 | 2019-03-22 | 北京恒泰实达科技股份有限公司 | 基于静止卫星监测资料的输电线路山火监测预警方法 |
CN109632106A (zh) * | 2019-01-03 | 2019-04-16 | 北京师范大学 | 一种遥感地表温度产品角度效应校正方法 |
CN110516816A (zh) * | 2019-08-30 | 2019-11-29 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于机器学习的全天候地表温度生成方法及装置 |
CN110866467A (zh) * | 2019-10-30 | 2020-03-06 | 核工业北京地质研究院 | 一种航空中红外高光谱数据温度和发射率反演方法 |
CN111563228A (zh) * | 2020-05-07 | 2020-08-21 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于地表入射短波辐射的山地地表反射率地形改正方法 |
-
2020
- 2020-10-26 CN CN202011153753.8A patent/CN112487346B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR101207925B1 (ko) * | 2010-02-22 | 2012-12-04 | 공주대학교 산학협력단 | 지표면 온도 산출 방법 |
WO2013077557A1 (ko) * | 2011-11-21 | 2013-05-30 | 한국항공우주연구원 | 성층권 비행선과 반사판을 이용한 지표면의 온도 조절 방법 |
CN103810387A (zh) * | 2014-02-13 | 2014-05-21 | 中国科学院地理科学与资源研究所 | 基于modis数据的地表蒸散发全遥感反演方法及系统 |
CN104657935A (zh) * | 2015-03-09 | 2015-05-27 | 广西壮族自治区气象减灾研究所 | 一种近地面气温反演方法 |
CN109509319A (zh) * | 2018-12-29 | 2019-03-22 | 北京恒泰实达科技股份有限公司 | 基于静止卫星监测资料的输电线路山火监测预警方法 |
CN109632106A (zh) * | 2019-01-03 | 2019-04-16 | 北京师范大学 | 一种遥感地表温度产品角度效应校正方法 |
CN110516816A (zh) * | 2019-08-30 | 2019-11-29 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于机器学习的全天候地表温度生成方法及装置 |
CN110866467A (zh) * | 2019-10-30 | 2020-03-06 | 核工业北京地质研究院 | 一种航空中红外高光谱数据温度和发射率反演方法 |
CN111563228A (zh) * | 2020-05-07 | 2020-08-21 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于地表入射短波辐射的山地地表反射率地形改正方法 |
Non-Patent Citations (1)
Title |
---|
赵伟等: "基于Landsat8热红外遥感数据的山地地表温度地形效应研究", 《遥感技术与应用》 * |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113269464A (zh) * | 2021-06-10 | 2021-08-17 | 中国科学院地理科学与资源研究所 | 一种生态恢复评估方法和生态恢复评估装置 |
CN113269464B (zh) * | 2021-06-10 | 2024-04-23 | 中国科学院地理科学与资源研究所 | 一种生态恢复评估方法和生态恢复评估装置 |
CN113449438A (zh) * | 2021-07-30 | 2021-09-28 | 北京环境特性研究所 | 一种可见光至热红外波段对地遥感图像仿真方法 |
CN113449438B (zh) * | 2021-07-30 | 2023-10-10 | 北京环境特性研究所 | 一种可见光至热红外波段对地遥感图像仿真方法 |
CN113743000A (zh) * | 2021-08-13 | 2021-12-03 | 电子科技大学 | 一种生成高时间分辨率全天候地表温度的方法 |
CN113743000B (zh) * | 2021-08-13 | 2023-04-18 | 电子科技大学 | 一种生成高时间分辨率全天候地表温度的方法 |
CN113887780A (zh) * | 2021-08-26 | 2022-01-04 | 国家卫星气象中心(国家空间天气监测预警中心) | 一种卫星遥感地表温度的估计方法、装置及设备 |
CN113887780B (zh) * | 2021-08-26 | 2023-11-24 | 国家卫星气象中心(国家空间天气监测预警中心) | 一种卫星遥感地表温度的估计方法、装置及设备 |
CN113945282A (zh) * | 2021-09-03 | 2022-01-18 | 中国空间技术研究院 | 红外遥感卫星温度反演精度指标分配系统及方法 |
CN114152350A (zh) * | 2021-12-09 | 2022-03-08 | 中国农业科学院农业资源与农业区划研究所 | 一种考虑城市三维几何结构影响的地表温度反演方法 |
CN114152350B (zh) * | 2021-12-09 | 2022-07-01 | 中国农业科学院农业资源与农业区划研究所 | 一种考虑城市三维几何结构影响的地表温度反演方法 |
CN115374390A (zh) * | 2022-10-26 | 2022-11-22 | 中国科学院地理科学与资源研究所 | 一种高光谱热红外地表温度和发射率波谱协同反演方法 |
Also Published As
Publication number | Publication date |
---|---|
CN112487346B (zh) | 2021-07-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112487346B (zh) | 一种山地地表温度遥感反演方法 | |
Tonooka | An atmospheric correction algorithm for thermal infrared multispectral data over land-a water-vapor scaling method | |
CN109883957B (zh) | 基于modis影像的表观反射率模型构建方法、系统及定标方法 | |
Gao et al. | Refinement of wavelength calibrations of hyperspectral imaging data using a spectrum-matching technique | |
CN107655574B (zh) | 时序热红外遥感地表温度/发射率同时反演的方法及系统 | |
CN102736128A (zh) | 无人机光学遥感影像数据处理方法及装置 | |
CN106778516B (zh) | 一种利用中国风云三号卫星遥感数据确定地表温度的方法及装置 | |
CN111323129A (zh) | 一种基于宽波段热红外影像的地表温度反演方法 | |
Zhou et al. | An improved temperature and emissivity separation algorithm for the advanced Himawari imager | |
CN109632106B (zh) | 一种遥感地表温度产品角度效应校正方法 | |
CN110866467A (zh) | 一种航空中红外高光谱数据温度和发射率反演方法 | |
CN114152350B (zh) | 一种考虑城市三维几何结构影响的地表温度反演方法 | |
Nie et al. | Land surface temperature and emissivity retrieval from nighttime middle-infrared and thermal-infrared Sentinel-3 images | |
Yu et al. | Scale mismatch between in situ and remote sensing observations of land surface temperature: Implications for the validation of remote sensing LST products | |
CN116519557A (zh) | 一种气溶胶光学厚度反演方法 | |
Zhang et al. | Development of the direct-estimation albedo algorithm for snow-free Landsat TM albedo retrievals using field flux measurements | |
Zeng et al. | Land surface temperature and emissivity retrieval from nighttime middle and thermal infrared images of Chinese Fengyun-3D MERSI-II | |
Tonooka et al. | Validation of ASTER/TIR standard atmospheric correction using water surfaces | |
CN114296061A (zh) | 基于多元变量检测和不同辐射传输模型的交叉定标方法 | |
Montes et al. | Comparing landsat-7 etm+ and aster imageries to estimate daily evapotranspiration within a mediterranean vineyard watershed | |
CN115269549A (zh) | 一种耦合物理-统计-深度学习的大气水汽反演方法 | |
Zhou et al. | Land surface albedo estimation with Chinese GF-1 WFV data in Northwest China | |
Le Borgne et al. | Validation of the OSI SAF radiative fluxes | |
Czapla-Myers et al. | Implication of spatial uniformity on vicarious calibration using automated test sites | |
CN114003852B (zh) | 一种地表下行太阳辐射通量的估算方法和装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |