CN103983360A - 一种基于hj-1b irs卫星数据反演lst的方法 - Google Patents

Abstract

本发明针对国产HJ-1B卫星IRS传感器数据提供一种地表温度(LST)反演技术，该技术主要包括以下步骤：利用MODIS的植被指数产品反演地表发射率(LSE)；利用单窗算法反演LST；利用遥感图像时空融合技术将不同时相的LST数据变换到基准时相下；区域数据镶嵌与裁切。本发明利用遥感图像时空融合技术，通过国外MODIS传感器高精度的LST产品来提高HJ-1B IRS数据LST反演的精度，同时实现了生产大区域同时相HJ-1B IRS LST产品的目的。

Description

一种基于HJ-1B IRS卫星数据反演LST的方法

技术领域

本发明涉及一种基于HJ-1B IRS卫星数据反演LST的方法，该方法是一种利用国产HJ-1B卫星IRS传感器数据反演地表温度LST的实用技术方法，能够应用在农业和环境监测等领域，属于卫星遥感技术领域。

背景技术

陆地表面温度(Land Surface Temperature，LST)是一个重要的地球物理参数。利用卫星遥感数据提取海洋温度(Sea Surface Temperature，SST)技术已较为成熟，可以在全球范围内达到1K的精度。由于陆地表面比海洋表面复杂得多，导致陆地表面温度反演的精度较低，陆地表面温度反演成了一个亟待解决的科学难题。

随着遥感技术的发展，对地热红外遥感在环境质量与灾害监测、气象气候热场监测、农业干旱监测、城市热岛调查等方面发挥着重要的作用，为了满足我国不断发展的热红外定量化应用需求，我国于2008年9月6日在太原卫星发射中心成功地发射了环境一号A、B卫星(HJ-1A、B卫星)。其中HJ-1B上搭载了红外相机IRS，热红外通道的光谱范围是1015～1215μm，星下点空间分辨率为300m，量化等级为10bit，幅宽为720km。但由于存在传感器本身温控技术水平低、绝对定标系数不精确等问题，导致HJ-1B IRS数据的LST反演精度不高、不稳定，严重制约了其定量化应用。

为了提高HJ-1B IRS数据LST反演精度，国内科学家已经做了很多有价值的工作，有学者以发射前实验室定标为基础，利用星上黑体定标数据对HJ-1B热红外通道进行星上定标，有学者以青海湖不同时相的美国TERRA卫星MODIS传感器31通道数据对HJ-1B热红外通道进行交叉定标，提高了HJ-1B热红外通道定标系数的精度，还有学者根据HJ-1B IRS数据的特点开发了专门的单窗LST反演算法，这些工作使得反演的LST数据更加精确，但距离国外传感器反演的LST数据精度还有一定差距，比如美国TERRA卫星MODIS传感器的LST数据产品。

国外科研人员利用某些地物在热红外波段比辐射率的稳定性和大气层顶太阳中波波段的辐射率变化极小的现象，利用MODIS的31和32两个热红外波段，开发了白天/黑夜适用的推广的分裂窗口LST算法，该算法通过实地测量验证，在晴空条件下精度可达到1K。由于TERRA卫星与HJ-1B卫星轨道设计相似，两者成像时间通常相差半个小时以内，可以假定两种数据反演的LST数值上将会非常接近，并且TERRA卫星重访周期短，同一地区每天都有成像数据，所以对于HJ-1B IRS数据几乎都可以找到对应的MODIS数据，因此本发明的出发点便是利用MODIS LST数据来提高HJ-1B IRS数据LST反演的精度。

发明内容

本发明的目的是针对HJ-1B IRS数据在定量遥感中的应用，提供一种有效的LST反演技术。本技术利用遥感图像时空融合技术，通过国外传感器MODIS高精度的LST反演产品来提高HJ-1B IRS数据LST反演的精度，同时实现了大区域内不同实相HJ-1B IRS数据反演LST到同一时像的目的，是一种针对国产HJ-1B IRS数据反演地表温度LST的实用技术方法。

本发明的基本思路为：对于单景或多景获取时间在一个月以内的、地理位置相近的多景HJ-1B IRS数据，利用单窗LST算法反演对应时像的地表温度，根据每景HJ-1B IRS数据的获取时间与地理范围下载对应的美国TERRA卫星MODIS传感器数据反演的LST产品MOD11A1/2数据，然后利用遥感图像时空融合的办法将不同时像的多景HJ-1B IRS LST数据逐像素变换到一个固定的时像下。本发明是一种利用HJ-1B IRS数据生产大区域的、固定时像下的LST产品的技术。

本发明是一种从HJ-1B IRS数据反演LST的方法，该方法具体步骤如下：

步骤一、获取时间在一个月以内的、地理位置相近的n景HJ-1B IRS数据，其中，n≥1，根据获取时间及地理位置下载对应的MODIS的植被指数产品，利用该产品反演每景HJ-1B IRS数据对应的地表发射率(Land Surface Emissivity，LSE)；

步骤二、利用单窗算法反演每景HJ-1B IRS数据的LST；

步骤三、下载与每景HJ-1B IRS数据获取时间及地理位置对应的MODIS地表温度产品，并选定某一时像作为基准时相，利用遥感图像时空融合技术将步骤二中获取的不同时相的HJ-1B IRS LST数据变换到基准时相下；

步骤四、将变换到基准时相下的n景HJ-1B IRS LST数据镶嵌成一幅图，并根据矢量数据裁切出该区域的LST反演结果。

步骤一中所述对于n(n≥1)景获取时间在一个月以内的、地理位置相近的多景HJ-1B IRS数据，处理的对象是多景数据主要是针对定量遥感对LST反演的具体需求而言的，实际应用中研究区域通常是大区域，比如真个中国区域，或者省级行政区域，通常需要多景HJ-1B IRS数据才能完成覆盖，为了保证多景HJ-1B IRS数据地表温度的相似性，限定数据的获取时间在一个月以内。步骤一中所述的MODIS的植被指数产品为与HJ-1B卫星成像时间接近的TERRA卫星获取的数据生产的，目前最适用的产品为16天合成的MOD13Q1植被指数产品，标识名称：MODIS/Terra Vegetation Indices16-Day L3Global250m S1N GridV005。步骤一中所述的反演每景HJ-1B IRS数据对应的地表发射率(Land Surface Emissivity，LSE)，反演的方法采用Van de Griend and Owe’s(1993)提出的NDVI域值法。

步骤二中所述单窗算法基于Jiménez-and Sobrin(2003)提出的针对Landsat TM/ETM+数据第6波段的LST反演算法，经过扩展应用于HJ-1B IRS数据的LST反演。

步骤三中所述的与每景HJ-1B IRS数据获取时间及地理位置对应的MODIS地表温度产品，优先选用与HJ-1B IRS数据同一天成像的，此时成像时间相差半小时以内，对应的数据名称为MODIS/TerraMOD11A1Daily daytime1km SIN Grid LST，由于该数据部分区域可能出现缺失，对于缺失的部分可以采用时间相近的8天合成LST产品，名称为MODIS/Terra MOD11A28-Day daytime1km SIN Grid LST。步骤三中所述的遥感图像时空融合技术基于Gao et al.(2006)提出的融合框架spatial and temporal adaptivereflectance fusion model(STARFM)，该框架原本针对Landsat TM/ETM+第1～5、7波段地表反射率数据与对应的MODIS数据的融合，经过扩展应用于HJ-1B IRS LST数据与对应的MODIS数据的融合。

步骤四是根据定量遥感应用的实际情况设计的，目的是为了获得某一大区域固定时像下的LST数据，并根据HJ-1B IRS数据的特性开发了包含融合的数据镶嵌算法。

本发明与现有技术相比有如下特点：克服了HJ-1B IRS由于传感器本身不够成熟导致的LST反演精度低且不稳定的问题，本发明依赖MODIS数据进行HJ-1B IRS数据的LST反演，由于MODIS传感器先进，且其采用的LST劈窗算法精度高，本发明HJ-1B IRS LST反演结果在稳定性与精度上都优于现有的单窗算法结果。本发明对多景数据采用遥感时空融合技术实现了生产大区域、同时像的LST产品，更加贴近实际目前定量遥感对大区域目标LST研究的需求，扩展了HJ-1B IRS数据LST反演的定量应用价值。

附图说明：

图1是HJ-1B IRS数据LST反演流程图

图2是地表发射率LSE与植被指数NDVI之间的函数关系示例图

图3是陕西省LST时空融合结果图

图中符号说明如下：

图1中，t₁，t₂…t_n表示多景数据不同的成像时间，t_i表示遥感图像时空融合的目标基准时刻。

具体实施方式：

采用本发明实现HJ-1B IRS卫星数据反演LST的实施例如图1所示，现结合附图对其进行描述。

图1所示反演流程的输入数据为同区域内多景HJ-1B IRS卫星数据，该数据IRS传感器有4个波段，可用于LST反演的是波长位于10.5-12.5μm的第四波段，空间分辨率300米，幅宽720千米，最大视场角30°。数据获取时间不同，为t₁，t₂，…t_n时刻，限定在1个月以内，根据HJ-1B IRS卫星数据下载对应的MODIS数据，具体为16天合成的MOD13Q1植被指数产品及每天或8天合成的MOD11A1/A2地表温度产品。MODIS数据通过网络下载，网址：http://reverb.echo.nasa.gov/reverb/。

图1所示地表发射率估算。地表发射率(LSE)即比辐射率，是物体的出射度与同温度黑体出射度之比。利用同时像MODIS的NDVI数据估算LSE，按照数值将NDVI数据分成6个分段，当NDVI值位于0.157到0.727时，根据Griend and Owe’s(1993)的研究结果，NDVI与LSE的关系曲线可以利用一个对数函数拟合，具体见图2，地表发射率ε由以下公式计算：

ε＝1.0094+0.047×ln(NDVI)

对于NDVI值低于0.157的象素，根据经验地表主要为裸地和水体，地表发射率取值0.955，对于NDVI超过0.727的象素，认为全部为植被覆盖，地表发射率取值0.99。

图1所示地表温度反演。HJ-1B IRS数据LST反演采用单窗算法，该单窗算法基于Jiménez-andSobrin(2003)提出的针对Landsat TM/ETM+数据第6波段的LST反演算法框架，并根据HJ-1B IRS第四波段数据特征简化计算步骤、拟合关键参数得到。记T_S为LST，单位为K，则该算法计算公式为：

T_S＝γ[ε^-1(ψ₁L_sensor+ψ₂)+ψ₃]+δ

其中：

$\mathrm{\γ}={\left\{\frac{c^2{L}_{\mathrm{sensor}}}{T_{\mathrm{sensor}}^2}\right[\frac{{\mathrm{\λ}}^4}{c_1}{L}_{\mathrm{sensor}}+{\mathrm{\λ}}^{-1}\left]\right\}}^{-1}$

δ＝-γL_sensor+T_sensor

其中地表发射率ε由图1所示地表发射率估算步骤获得，L_sensor为感器接收到的辐射亮度，单位W*m^-2*sr^-1*μm^-1，该值根据辐射定标公式，由HJ-1B IRS原始数据的DN值(digital number)计算得到，公式如下：

L_sensor＝DN/gain+bias

其中偏移和增益参数，HJ-1B IRS第四波段的参数gain＝61.472，bias＝-44.598。

T_sensor为对应的大气层顶亮度温度，由L_sensor计算出，公式如下：

$T_{\mathrm{sensor}}=\frac{c_1{\mathrm{\λ}}_c^{-1}}{\ln (c_2{\mathrm{\λ}}_c^{-5}/L_{\mathrm{sensor}}+1)}$

其中c₁＝1.43876869×10⁴与c₂＝1.1909561×10⁴，λ_c计算公式如下：

${\mathrm{\λ}}_c=\frac{{\∫}_{\mathrm{\λ}\min }^{\mathrm{\λ}\max }f\left(\mathrm{\λ}\right)\mathrm{\λd\λ}}{{\∫}_{\mathrm{\λ}\min }^{\mathrm{\λ}\max }f\left(\mathrm{\λ}\right)\mathrm{d\λ}}$

其中f(λ)为HJ-1B IRS第4波段，即热红外波段的光谱响应函数。

ψ₁、ψ₂、ψ₃为大气参数，由大气水汽含量ω计算得到，采用Duan et al(2008)针对HJ-1B IRS数据拟合出的结果，公式如下：

ψ₁＝0.0412ω²+0.0936ω+0.9856

ψ₂＝-0.7174ω²-0.8812ω-0.3941

ψ₃＝0.2639ω²+0.6499ω+0.4703

大气水汽含量ω的值根据HJ-1B IRS数据获取的季节、地区不同而定，是LST反演的关键参数，传统获取ω值的方法有两种：一是利用遥感数据相应波段反演，比如MODIS利用波段在0.86-1.20um的第17、18、19波段反演水汽含量；二是根据研究区域情况人工给定一个经验值。第二种方式虽然存在相对大的不确定性，但在实际研究中较为常用，本技术采用人工给定经验值的方式，虽然会引入LST反演的不确定性，但对最终结果精度的影响不大，因为在图1所示遥感数据时空融合步骤中，LST反演的最终精度主要由MODIS LST产品决定。

图1所示遥感图像时空融合。该步骤目的是将不同时相t₁，t₂，…t_n的HJ-1B IRS数据统一到特定时相t_i下，达到获取大区域内多景HJ-1B IRS数据同一时刻LST数据的目的，该步骤时空融合的参考数据为与每景HJ-1B IRS数据对应的MODIS LST产品数据，采用的融合模型基于Gao et al.(2006)提出的STARFM模型(spatial and temporal adaptive reflectance fusion model，STARFM)，在STARFM的基础上通过步骤简化获得。该模型公式如下：

$F(x_{w/2},y_{w/2},t_2)=\underset{i=1}{\overset{w}{\mathrm{\Σ}}}\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{W}_{\mathrm{ijk}}\×(C(x_i,y_j,t_2)+F(x_i,y_j,t_k)-C(x_i,y_j,t_k))$

其中F表示高分辨率遥感图像中的象素点，这里指HJ-1B IRS数据，C表示低分辨率遥感图像中的象素点，这里指MODIS LST产品数据，(x_i，y_j)表示象素位置，t_k，k＝1，2，…n表示遥感图像不同的获取时间，t₂表示需要归一化到的时间，w为移动窗口的大小，(x_w/2，y_w/2)表示移动窗口中的中心象素点。W_ijk为该时空融合模型中的关键权重函数，计算公式如下：

$W_{\mathrm{ijk}}=\frac{1}{E_{\mathrm{ijk}}}/\underset{i=1}{\overset{w}{\mathrm{\Σ}}}\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(1/E_{\mathrm{ijk}})$

其中参数E_ijk的计算公式如下：

E_ijk＝ln(S_ijk×A+1)×ln(T_ijk×B+1)×D_ijk

其中A与B为尺度因子，根据经验这里取值1000，S_ijk表示象素的光谱差异，即象素值的差异。T_ijk表示象素值的时相差异，即获取时间的差异。D_ijk表示象素位置的几何距离的差异。可见参数E_ijk的值由移动窗口w所有象素的象素值、获取时间以及距离中心象素点的距离决定的。HJ-1B IRS数据的空间分辨率为300米，MODIS LST产品空间分辨率1000米，为了实现数据的象素点对应，采用三次卷积插值的方式将MODIS LST产品重采样到空间分辨率300米。

图1所示区域数据镶嵌、裁切。该步骤对于生产多景HJ-1B IRS LST覆盖的大区域同时相产品是必要的，将多景经过遥感图像时空融合的LST结果数据按照统一的地图投影镶嵌成一幅图像，并根据生产区域的矢量文件将对应区域的数据裁切出来。该步骤利用现有大多数商用软件即可完成，比如ENVI、ERDAS等都提供相应的镶嵌与数据裁切工具。传统的镶嵌软件对于相邻两景数据重叠区域的处理只涉及数据的取舍问题，而不考虑重叠区域数据的融合，这对于具有固定重叠区域的遥感数据是适用的，裁切过程中丢弃的重叠区域占总体数据的比例通常不大，但对于HJ-1B IRS数据就可能会损失很大比例的数据信息，因为HJ-1B卫星是非回归轨道，数据之间的重叠区域大小不固定，当两景数据重叠区域超过图像的一半时，直接镶嵌会丢弃大量数据信息，尤其是当重叠区域存在云与云下阴影时，针对这个问题，本发明设计了一个考虑数据融合的HJ-1B IRS数据镶嵌算法。该算法假定在一个区域内有m景数据重叠，数据获取的时间有n个，分别为t₁，t₂，…t_n，m景数据在进行遥感图像时空融合处理的目标时间为t_i。记P_k∈m(x，y，t_j)表示重叠区域m景数据中的第k景数据，数据获取时间为t_j，按照地理坐标对应的象素位置(x，y)的象素值，则镶嵌融合后的象素值P(x，y，t_i)的计算公式为：

$P(x,y,t_i)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k{P}_k(x,y,t_j)$

其中α_k，k∈[1，m]为权重系数，且有α_k的计算过程如下：

首先计算m个象素点P_k∈m(x，y，t_j)的相似性。方法是以象素位置(x，y)为中心计算w×w象素尺寸的窗口，w可以取值8或者16等值。然后统计窗口象素属于重叠区域内并且不是背景值的象素，对这些象素位置计算m景图像两两窗口之间的相关系数，相关系数低于一定域值表明当前象素位置(x，y)在其中一景图像中可能为云或云下阴影，或者地表发生变化，若m≥3，则设定与其它数据相关系数都低的那景数据对应的权重α_k＝0，若m＝2，则设定与数据获取时间距离基准时间t_i远的那景数据对应的权重α_k＝0。对于相关性高的象素点，α_k权重分值根据数据获取时间确定，距离基准时间t_i约近α_k值越大，且有 $\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k=1.$

图1所示整个处理流程的一个实例已经在PC平台上实现，支持Windows与Linux操作系统，目前支持大多数HJ-1B IRS LST数据的反演，具有生产大区域同时相、空间分辨率300米的LST产品的能力，本发明计算原理简单，计算量小，处理速度快，结合简单的OpenMP加速技术后，在普通单机32位、Intelcore i7-3770，CPU3.4GHz配置下，对于单景HJ-1B IRS数据平均处理时间小于5分钟，软件鲁棒性较强，满足HJ-1B IRS LST数据生产自动、快速处理的工程化需求。

图3所示陕西省LST时空融合结果图，该图是本发明实际生产大区域同时相LST产品的一个具体实例，利用6景HJ-1B IRS数据，数据的获取时间有4个，利用对应的MODIS LST产品，得到了陕西省统一时相的LST产品。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。

Claims (5)

1.一种基于HJ-1B IRS卫星数据反演陆地表面温度的方法，该方法是一种针对国产HJ-1B IRS数据反演陆地表面温度(Land Surface Temperature，LST)的实用技术方法，其特征在于包含以下步骤：

步骤一、获取时间在一个月以内的、地理位置相近的n景HJ-1B IRS数据，其中n≥1，根据获取时间及地理位置下载对应的MODIS的植被指数产品，利用该产品反演每景HJ-1B IRS数据对应的地表发射率(Land Surface Emissivity，LSE)；

步骤二、利用单窗算法反演每景HJ-1B IRS数据的LST；

步骤三、下载与每景HJ-1B IRS数据获取时间及地理位置对应的MODIS地表温度产品，并选定某一时相作为基准时相，利用遥感图像时空融合技术将步骤二中获取的不同时相的HJ-1B IRS LST数据变换到基准时相下；

步骤四、将变换到基准时相下的n景HJ-1B IRS LST数据镶嵌成一幅图，并根据矢量数据裁切出该区域的LST反演结果。

2.根据权利要求1中所述的方法，其特征在于：步骤一中所述的反演每景HJ-1B IRS数据对应的地表发射率(Land Surface Emissivity，LSE)的方法采用NDVI域值法。

3.根据权利要求1中所述的方法，其特征在于：步骤二中所述的单窗算法为在适用于Landsat TM/ETM+数据第6波段的LST反演算法的基础上，经过扩展得到的HJ-1B IRS数据的LST反演算法。

4.根据权利要求1中所述的方法，其特征在于：步骤三中所述的遥感图像时空融合技术采用STARFM融合框架(spatial and temporal adaptive reflectance fusion model)，HJ-1B IRS数据的空间分辨率为300米，MODIS LST产品空间分辨率1000米，为了实现数据的象素点对应，采用三次卷积插值的方式将MODISLST产品重采样到空间分辨率300米。

5.根据权利要求1中所述的方法，其特征在于：步骤三中所述的变换到基准时相下的HJ-1B IRS LST数据镶嵌成一幅图的方法采用的是一种考虑数据融合的HJ-1B IRS数据镶嵌算法，该算法假定在一个区域内有m景数据重叠，数据获取的时间有n个，分别为t₁，t₂，…t_n，m景数据在进行遥感图像时空融合处理的目标时间为t_i，记P_k∈m(x，y，t_j)表示重叠区域m景数据中的第k景数据，数据获取时间为t_j，按照地理坐标对应的象素位置(x，y)的象素值，则镶嵌融合后的象素值P(x，y，t_i)的计算公式为：

$P(x,y,t_i)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k{P}_k(x,y,t_j),$

其中α_k，k∈[1，m]为权重系数，且有α_k的计算过程如下：

首先计算m个象素点P_k∈m(x，y，t_j)的相似性，方法是以象素位置(x，y)为中心计算w×w象素尺寸的窗口，w可以取值8或者16等值。然后统计窗口象素属于重叠区域内并且不是背景值的象素，对这些象素位置计算m景图像两两窗口之间的相关系数，相关系数低于设定域值，则表明当前象素位置(x，y)在其中一景图像中为云、云下阴影或者地表发生变化，若m≥3，则设定比其它数据相关系数都低的那景数据对应的权重α_k＝0，若m＝2，则设定与数据获取时间距离基准时间t_i远的那景数据对应的权重α_k＝0，对于相关性高的象素点，α_k权重分值根据数据获取时间确定。