一种基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法

1.本发明涉及参数遥感反演和背景的红外辐射特性研究领域，具体涉及一种基于热红外遥感反演雪表温度和基于matlab程序计算雪表红外辐射量的方法。


背景技术：

2.冰冻圈通常认为由地球表层每年至少部分时间处在零度以下所形成的各类积雪、冰川、河流湖泊中的淡水冰、海冰、地下冰及多年冻土等部分组成。雪在冰冻圈中分布最为广阔，是地球表面最为活跃的自然要素之一，也是全球变化长期监测的自然资源与环境要素之一。积雪的特性(如积雪面积、雪深、雪表温度、积雪反射率等)是全球水能平衡模型中的主要输入参数。对于大陆尺度而言，大范围雪盖可影响气候变化、地表辐射平衡与能量交换、水循环等；而对于局地和流域尺度而言，雪盖可影响天气、工农业和生活用水资源、环境、寒区工程等一系列与人类活动有关的要素。
3.由于受到纬度、海拔、植被覆盖、冻土分布和土壤含水量等时空多变因素的影响，积雪表层温度的时空分布模式很复杂，使得获取宏观范围的雪温资料很困难。目前雪温资料主要来自于气象站点的观测和人工测量，但是站点观测亦或是人工测量只能提供空间上离散的有限点观测数据，而地球系统模型往往需要空间上连续的面数据而非点数据。因此，传统的气象资料需要通过如反距离权重插值、克里金插值等空间插值的再处理方法得到所需的面数据，但是在站点密度较低的地区，该再处理方法得到的数据往往满足不了精度要求。随着全球变化研究的不断深入，常规气象站点观测以及人工测量已经无法满足地球系统科学研究的需要。雪温时空分布资料的缺乏在一定程度上制约了地学模型的发展，从而影响我们对全球气候和环境的认识。与积雪表层温度息息相关的雪表面的红外辐射通量，在大范围内对其进行准确的计算和场景展现，对雪地背景中目标的隐身与识别具有非常重要的指导意义。
4.近些年来，遥感技术的飞速发展为快速地获取大尺度的雪温时空分布信息提供了新的途径和方法。卫星遥感最突出的优势在于能够提供大范围的空间和时间上连续的雪表及大气信息，基于遥感影像反演得到的雪温数据提供了比站点观测数据更理想的空间异质度信息。此外，遥感还可以根据历史存档数据提取过去时相的气温，着对于过去气候、生态、植被和冰川冻融等研究具有重要的意义。
5.利用热红外遥感反演地表温度的方法已经比较成熟，而积雪表层温度的反演研究开展得很少，目前应用最普遍的雪表层温度遥感反演方法是多通道算法，该方法最初被应用于地温遥感反演上，后通过在低温范围内对plank函数进行了线性化，修正模型参数使其能够应用于积雪表温的反演研究上。在积雪表层温度的反演研究中，没有给出涉及输入参数具体的估算方案，也没有针对性地对高纬度高海拔的常年积雪进行研究；在积雪表面的红外特性研究中，也没有合适的计算模型可以将热红外遥感反演的结果计算为大场景范围内雪表面的红外辐射通量。


技术实现要素：

6.本发明目的是针对现有积雪表温遥感反演技术存在的不足和雪表面红外辐射计算模型的空缺，在考虑地表能量平衡过程的基础上，综合选用了积雪指数ndsi、地表反照率、大气水汽含量、大气透过率、地表比辐射率、亮度温度等多个地表与气候因子来建立其与积雪表层温度之间的关系，并确立了各输入参数的估算方案，建立了积雪表温的遥感反演模型；基于积雪表温反演模型的反演结果，利用matlab程序构建了积雪表面红外辐射通量的计算模型，完成了大场景积雪表面红外特性的研究。
7.为实现上述目的，本发明提供如下技术方案：
8.本发明提出一种基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法包括如下步骤：
9.(1)建立了积雪表面温度的遥感反演模型：
10.ts＝b0+b1t
31
+b2t
32 (1)
11.其中：
[0012][0013]c31
＝τ
31
ε
31
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀd31
＝(1-τ
31
)[1+τ
31
(1-ε
31
)]
[0014]c32
＝τ
32
ε
32
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀd32
＝(1-τ
32
)[1+τ
32
(1-ε
32
)]
[0015]
式中，τ
31
和τ
32
分别是31波段和32波段的大气透过率，ε
31
和ε
32
分别是31波段和32波段的地表比辐射率，t
31
和t
32
分别是31波段和32波段对应影像的亮度温度，ts为积雪表面的温度。
[0016]
(2)确立了雪表温遥感反演模型各输入参数的估算方案：
[0017]
首先是大气透过率的估算方案：
[0018]
在中纬度地区夏季时，τ
31
＝-0.10671ω+1.04015r2＝0.995
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0019]
τ
32
＝-0.12577ω+0.99229 r2＝0.996 (3)
[0020]
在中纬度地区冬季时，τ
31
＝-0.1041ω+0.92314r2＝0.996
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0021]
τ
32
＝-0.13722ω+0.97686 r2＝0.987 (5)
[0022]
式中，ω＝{[α-ln(ρ
19
/ρ2)]/β}2ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0023]
ρ
19
和ρ2分别是modis影像的第19和2波段的反射率，ω为大气水汽含量，单位为g/cm2。对于混合型地表α＝0.02，β＝0.651；裸土地表α＝-0.04，β＝0.651；植被覆盖地表α＝0.0012，β＝0.651。
[0024]
其次是地表比辐射率的估算方案：
[0025]
该估算方法建立在modis雪盖率反演图基础上，地表覆盖物的光谱反射率均来自于aster光谱库。基本思想是首先是对像元内的积雪覆盖率进行判断，辅助以mcd12q1土地覆盖类型数据进行分析，然后对应不同的地表组成情况赋予不同的比辐射率值。当像元的积雪覆盖度为0时，认为该像元无积雪覆盖，则令其地表比辐射率为0；当像元的积雪覆盖度为1时，认为该像元时纯雪像元，则令其地表比辐射率为aster光谱库积雪光谱数据的积分计算值；当像元的积雪覆盖度为0～1时，根据mcd12q1遥感资料，若像元为林区，认为像元是积雪和植被的混合像元，若像元为非林区，认为像元是积雪和裸土的混合像元。此时像元的
地表比辐射率可以利用雪盖率、积雪比辐射率、植被比辐射率和裸土比辐射率按比例计算得到。具体的计算公式如所示：
[0026][0027]
式中，ε
snow
、εv、εs分别是积雪、植被和裸土的比辐射率，可以利用aster光谱库中对应的光谱曲线按照相应波段积分计算得到。fsc为modis像元的积雪覆盖率，公式中的fsc
·
ε
snow
、(1-fsc)
·
εv和(1-fsc)
·
εs是积雪、植被和裸土在像元中的有效比辐射率，即为各端元比辐射率的面积加权值。
[0028]
最后是影像亮度温度的估算方案：
[0029][0030]
式中，li为地表的亮度温度；h为普朗克常量，取值6.626
×
10-34j·
s；c是光速取值为2.992
×
108m/s；k为玻尔兹曼常数，取值为1.3806
×
10-23
j/k；λ为波长，单位为m。
[0031]
(3)基于积雪表层温度的遥感反演模型和matlab计算程序构建积雪表面红外辐射通量的计算模型：
[0032]
在得到积雪的表面温度后，积雪表面的本身辐射通量可以从普朗克公式出发，通过对红外波段范围(对应大气窗口红外波段范围)积分得到：
[0033][0034]
式中，t为单元表面温度(k)；ε
snow
为积雪表面的比辐射率；λ2、λ1分别为红外波段范围的上下限(μm)；c1为第一辐射常数，等于3.742
×
108w
·
μm4/m2；c2为第二辐射常数，等于1.439
×
104μm
·
k。
[0035]
故雪地表面自身的辐射通量：
[0036][0037]
相对于现有技术：本发明提供了一种积雪表温的热红外遥感估算方法和积雪表面的红外辐射计算模型，能够快速准确的计算出积雪表面的温度和红外辐射通量的空间分布图。由于本方法是在考虑地表能量平衡过程的基础上通过回归分析方法建立的，因此同时具有理论模型和统计模型的优点，精度较高，且具有较强的实用性；
[0038]
本发明针对中国的西藏稳定积雪区进行了积雪表温的遥感反演，并与gldas雪温数据集进行对比分析，验证了反演方法的精度和可靠性；基于积雪表温的遥感反演结果，利用matlab程序构建了积雪表面的红外辐射计算模型。
附图说明
[0039]
图1：积雪表温观测值与公式(1)遥感估算值的折线图。
[0040]
图2：积雪表温冬季的遥感反演值(大场景)。
[0041]
图3：积雪表面冬季红外辐射通量计算值(3-5μm)。
[0042]
图4：积雪表面冬季红外辐射通量计算值(8-14μm)。
[0043]
图5：积雪表面红外辐射通量的计算模型的matlab程序。
[0044]
图6：热红外遥感反演雪表温度和计算雪表红外辐射量的流程图。
具体实施方式
[0045]
下面结合图6对本发明进一步详述：
[0046]
1)：建立积雪表面温度的遥感反演模型
[0047]
在考虑地表能量平衡过程的基础上，对地表热辐射传输方程中的大气上下行热辐射进行近似求解，并且利用大气平均作用温度代替大气向下平均作用温度，在低温范围内对planck函数进行线性化展开。综合选用了积雪指数ndsi、地表反照率、大气水汽含量、大气透过率、地表比辐射率、亮度温度等多个地表与气候因子来建立其与积雪表层温度之间的关系，并确立了各输入参数的估算方案，建立了积雪表温的遥感反演模型。
[0048]
ts＝b0+b1t
31
+b2t
32 (1)
[0049]
其中：
[0050][0051]c31
＝τ
31
ε
31
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀd31
＝(1-τ
31
)[1+τ
31
(1-ε
31
)]
[0052]c32
＝τ
32
ε
32
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀd32
＝(1-τ
32
)[1+τ
32
(1-ε
32
)]
[0053]
式中，τ
31
和τ
32
分别是31波段和32波段的大气透过率，ε
31
和ε
32
分别是31波段和32波段的地表比辐射率，t
31
和t
32
分别是31波段和32波段对应影像的亮度温度，ts为积雪表面的温度。
[0054]
2)：确立雪表温遥感反演模型各输入参数的估算方案
[0055]
通过modis遥感数据可以获取地表反射率、积雪指数、亮度温度等变量值，利用这些变量通过参数估算方案进一步计算到雪表温遥感反演模型需要的驱动参数，过程如下：
[0056]
a.亮度温度：首先将mod021km遥感数据通过在envi中的预处理得到研究区地表在31和32波段的辐亮度影像，再代入式(8)算出31和32波段的亮度温度。
[0057]
b.像元平均比辐射率：首先利用aster光谱库中积雪、植被和裸土的光谱数据，根据modis第31和32通道的宽度，采用积分方式计算得到中等雪粒子的比辐射率分别为0.98916和0.97983；植被的比辐射率为0.9851和0.9844；裸土的比辐射率为0.9832和0.9731。在此基础上，通过式(7)可以计算出像元平均比辐射率。
[0058]
c.大气透过率：首先利用mod021km遥感数据中的第2和19通道的反射率，通过两通道比值法计算出大气水汽含量；再利用中纬度冬季的大气透过率估算方程，计算大气透过率值。
[0059]
d.反演过程中像元中端元面积的比例依靠ndsi雪盖图来确定：若是积雪像元，其值为像元的积雪覆盖率fsc；若是植被或者裸土端元，则值为1-fsc。
[0060]
利用gldas雪温数据作为验证样本对遥感反演模型公式(1)进行精度验证，在长达一周的验证时间内误差稳定在2-4k范围内，表明本发明建立的遥感反演模型具有比较理想的精度。图1给出了遥感反演雪表温度值与雪温产品观测值的变化折线图。
[0061]
3)：构建积雪表面红外辐射通量的计算模型
[0062]
在得到积雪的表面温度后，积雪表面的本身辐射通量可以从普朗克公式出发，通过对红外波段范围(对应大气窗口红外波段范围)积分得到：
[0063][0064]
式中，t为单元表面温度(k)；ε
snow
为积雪表面的比辐射率；λ2、λ1分别为红外波段范围的上下限(μm)；c1为第一辐射常数，等于3.742
×
108w
·
μm4/m2；c2为第二辐射常数，等于1.439
×
104μm
·
k。
[0065]
故雪地表面自身的辐射通量：
[0066][0067]
考虑到雪表温遥感反演模型输出的是大场景的雪表温栅格数据，编写matlab程序，首先完成文件格式由“tif-ascii文件-tif”的转化，再构建出积雪表面在大气窗口3-5μm和8-14μm波段的红外辐射通量计算模型。
[0068]
选用2020年该年的modis l1b遥感数据，基于积雪表温的热红外遥感反演模型，完成四季的雪表温计算，并以栅格数据文件的格式，呈现大场景范围内雪表温的空间分布情况。图2给出了积雪表温四季的遥感反演值(大范围)。
[0069]
考虑大气窗口在红外波段范围有3-5μm和8-14μm两个波段，在得到积雪表温四季的遥感反演值基础上，利用编写的matlab程序分别计算出两个波段窗口积雪表面的红外辐射通量值。图3给出了积雪表面红外辐射通量计算值(3-5μm)；图4给出了积雪表面红外辐射通量计算值(8-14μm)；图5给出了积雪表面红外辐射通量的计算模型的matlab程序。
[0070]
该方法可以为研究区的雪表的温度遥感反演和雪表的红外辐射通量计算提供方法参考和方案支持，其中的积雪表面温度的遥感反演模型和各输入参数的估算方案可以基于modis等遥感数据快捷且准确地完成雪表的温度反演，为研究区的积雪温度监测提供数据支撑；其中的积雪表面红外辐射通量的计算模型可以较好地完成大气窗口下雪表的红外特性计算，并以栅格图像的形式输出结果，可以为雪地背景中目标的红外隐身和探测提供理论基础和方法支持。

技术特征：
1.一种基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，其特征在于，包括如下步骤：(1)确立了积雪表层温度遥感反演模型的基本形式，(2)建立了积雪表层温度的遥感反演模型和各驱动参数的确定方案，(3)基于积雪表层温度的遥感反演模型和matlab计算程序构建积雪表面红外辐射通量的计算模型。2.根据权利要求1所述的基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，其特征在于，建立了积雪表面温度的遥感反演模型具体步骤为：t
s
＝b0+b1t
31
+b2t
32
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)其中：其中：c
31
＝τ
31
ε
31 d
31
＝(1-τ
31
)[1+τ
31
(1-ε
31
)]c
32
＝τ
32
ε
32 d
32
＝(1-τ
32
)[1+τ
32
(1-ε
32
)]式中，τ
31
和τ
32
分别是31波段和32波段的大气透过率，ε
31
和ε
32
分别是31波段和32波段的地表比辐射率，t
31
和t
32
分别是31波段和32波段对应影像的亮度温度，t
s
为积雪表面的温度。3.根据权利要求1所述的基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，其特征在于，步骤2中，确立了雪表温遥感反演模型中大气透过率参数的估算方案：在中纬度地区夏季时，τ
31
＝-0.10671ω+1.04015 r2＝0.995
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)τ
32
＝-0.12577ω+0.99229 r2＝0.996
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)在中纬度地区冬季时，τ
31
＝-0.1041ω+0.92314 r2＝0.996
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)τ
32
＝-0.13722ω+0.97686 r2＝0.987
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)式中，ω＝{[α-ln(ρ
19
/ρ2)]/β}2ꢀꢀꢀꢀꢀ
(6)ρ
19
和ρ2分别是modis影像的第19和2波段的反射率，ω为大气水汽含量，单位为g/cm2；对于混合型地表α＝0.02，β＝0.651；裸土地表α＝-0.04，β＝0.651；植被覆盖地表α＝0.0012，β＝0.651。4.根据权利要求1所述的基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，其特征在于，步骤2中，确立了雪表温遥感反演模型中地表比辐射率参数的估算方案：式中，ε
snow
、εv、ε
s
分别是积雪、植被和裸土的比辐射率，可以利用aster光谱库中对应的光谱曲线按照相应波段积分计算得到；fsc为modis像元的积雪覆盖率，公式中的fsc
·
ε
snow
、(1-fsc)
·
εv和(1-fsc)
·
ε
s
是积雪、植被和裸土在像元中的有效比辐射率，即为各端元比辐射率的面积加权值。5.根据权利要求1所述的基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，其特征在于，步骤2中，确立了雪表温遥感反演模型中亮度温度参数的估算方案：
式中，l
i
为地表的亮度温度；h为普朗克常量，取值6.626
×
10-34
j
·
s；c是光速取值为2.992
×
108m/s；k为玻尔兹曼常数，取值为1.3806
×
10-23
j/k；λ为波长，单位为m。6.根据权利要求1所述的基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，其特征在于，基于积雪表层温度的遥感反演模型和matlab计算程序构建积雪表面红外辐射通量的计算模型：在得到积雪的表面温度后，积雪表面的本身辐射通量可以从普朗克公式出发，通过对红外波段范围(对应大气窗口红外波段范围)积分得到：式中，t为单元表面温度(k)；ε
snow
为积雪表面的比辐射率；λ2、λ1分别为红外波段范围的上下限(μm)；c1为第一辐射常数，等于3.742
×
108w
·
μm4/m2；c2为第二辐射常数，等于1.439
×
104μm
·
k；故雪地表面自身的辐射通量：7.根据权利要求6所述的基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，其特征在于，基于matlab程序完成计算结果从“tif文件-ascii文件-tif文件”的格式转化，初始的tif文件为预处理过的遥感产品数据结果，再转化为ascii文件的目的是便于读取像元单位的温度值并完成红外辐射通量的计算，再转化为tif文件的目的是将红外辐射通量的计算结果以栅格数据呈现，可视化的展现大场景积雪的红外辐射特性。

技术总结
本发明公布了一种基于热红外遥感反演雪表温度和计算雪表红外辐射量的方法，包括如下步骤：(1)确立了积雪表层温度遥感反演模型的基本形式，(2)建立了积雪表层温度的遥感反演模型和各驱动参数的确定方案，(3)基于积雪表层温度的遥感反演模型和Matlab计算程序构建积雪表面红外辐射通量的计算模型。本发明能够快速准确的计算出积雪表面温度和红外辐射量的空间分布图。由于本模型是在考虑地表能量平衡过程的基础上通过回归分析方法建立的，因此同时具有理论模型和统计模型的优点，精度较高，且具有较强的实用性。且具有较强的实用性。且具有较强的实用性。


技术研发人员：韩玉阁 陈兆炜
受保护的技术使用者：南京理工大学
技术研发日：2023.05.16
技术公布日：2023/8/16