CN110348107B - 一种云下像元真实地表温度重建方法 - Google Patents
一种云下像元真实地表温度重建方法 Download PDFInfo
- Publication number
- CN110348107B CN110348107B CN201910609171.7A CN201910609171A CN110348107B CN 110348107 B CN110348107 B CN 110348107B CN 201910609171 A CN201910609171 A CN 201910609171A CN 110348107 B CN110348107 B CN 110348107B
- Authority
- CN
- China
- Prior art keywords
- pixel
- cloud
- modis
- data
- lst
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 47
- 238000012545 processing Methods 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 26
- 239000002689 soil Substances 0.000 claims description 19
- 238000002310 reflectometry Methods 0.000 claims description 15
- 230000005855 radiation Effects 0.000 claims description 13
- 238000012937 correction Methods 0.000 claims description 9
- 238000003491 array Methods 0.000 claims description 4
- 125000004432 carbon atom Chemical group C* 0.000 claims description 4
- 239000000203 mixture Substances 0.000 claims 3
- 230000002354 daily effect Effects 0.000 description 47
- 238000012544 monitoring process Methods 0.000 description 12
- 238000004422 calculation algorithm Methods 0.000 description 11
- 238000011160 research Methods 0.000 description 10
- 238000010586 diagram Methods 0.000 description 9
- 230000008569 process Effects 0.000 description 9
- 230000007613 environmental effect Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 5
- 230000003595 spectral effect Effects 0.000 description 4
- 238000003908 quality control method Methods 0.000 description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 230000003203 everyday effect Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- CBENFWSGALASAD-UHFFFAOYSA-N Ozone Chemical compound [O-][O+]=O CBENFWSGALASAD-UHFFFAOYSA-N 0.000 description 1
- 239000000443 aerosol Substances 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 229910002056 binary alloy Inorganic materials 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 239000008264 cloud Substances 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 238000013524 data verification Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000004134 energy conservation Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 239000003595 mist Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/22—Matching criteria, e.g. proximity measures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computer Hardware Design (AREA)
- Geometry (AREA)
- Astronomy & Astrophysics (AREA)
- Remote Sensing (AREA)
- Multimedia (AREA)
- Image Processing (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
一种云下像元真实地表温度重建方法，其通过对连续多日的MODIS日产品数据进行处理，获得指定日期的MODIS日产品数据中的云下像元的真实地表温度，本发明所提供的一种云下像元真实地表温度重建方法，不仅能实现大面积云下像元LST的重建，且重建后的LST与地面实测的LST间具有较高的相关性。
Description
技术领域
本发明涉及遥感技术领域，特别涉及一种通过对MODIS数据进行处理，从而重建获得云下像元的真实地表温度的方法。
背景技术
地表温度(LST)常被定义为地球的表层温度。作为环境研究和资源管理的一个重要参数，被广泛应用于干旱监测、蒸散发、土壤水分估算和森林火灾探测等领域。20世纪70年代，遥感技术开始应用于地表温度反演中，其宏观、快速、经济等优势，弥补了传统地面监测在空间分布上的不足。高分辨辐射计(Advanced Very High Resolution Radiometer,AVHHR)和中分辨率成像光谱仪(Moderate Resolution Imaging Spectrometer,MODIS)是目前使用较为广泛的两种热红外传感器，在空间分辨率为1km的情况下，两者均能实现每日获取至少一次对地全覆盖观测。
1999年2月18日，美国成功地发射了地球观测系统(EOS)的第一颗先进的极地轨道环境遥感卫星Terra。它的主要目标是实现从单系列极轨空间平台上对太阳辐射、大气、海洋和陆地进行综合观测，获取有关海洋、陆地、冰雪圈和太阳动力系统等信息，进行土地利用和土地覆盖研究、气候季节和年际变化研究、自然灾害监测和分析研究、长期气候变率的变化以及大气臭氧变化研究等，进而实现对大气和地球环境变化的长期观测和研究的总体(战略)目标。2002年5月4日成功发射Aqua星后，每天可以接收两颗星的资料。搭载在Terra和Aqua两颗卫星上的中分辨率成像光谱仪(MODIS)是美国地球观测系统(EOS)计划中用于观测全球生物和物理过程的重要仪器。它具有36个中等分辨率水平(0.25um～1um)的光谱波段，每1-2天对地球表面观测一次。获取陆地和海洋温度、初级生产率、陆地表面覆盖、云、汽溶胶、水汽和火情等目标的图像。
MODIS数据涉及波段范围广，有36个离散光谱波段，光谱范围宽，从0.4μm到14.4μm全光谱覆盖。TERRA和AQUA卫星都是太阳同步极轨卫星，TERRA在地方时上午过境，AQUA在地方时下午过境。TERRA与AQUA上的MODIS数据在时间更新频率上相配合，加上晚间过境数据，对于接收MODIS数据来说可以得到每天最少2次白天和2次黑夜更新数据。这样的数据更新频率，对实时地球观测和应急处理(例如森林和草原火灾监测和救灾)有很大的实用价值。
MODIS仪器的地面分辨率包括250m、500m和1000m，扫描宽度为2330km。
MODIS陆地标准数据产品具有不同的时间分辨率，主要包括：每日产品(daily)、8天合成产品(8-Day)、16天合成产品(16-Day)、月合成产品(Monthly)、季度产品(Quarterly)以及年产品(Yearly)。
在过去几十年里，MODIS传感器已经获得了大量的地表热红外数据。经验证，基于普适性分裂窗算法和白昼算法计算反演的MODIS LST产品数据的误差小于1K，作为一种高质量的产品数据，MODIS LST产品数据被广泛应用到许多领域，尤其是干旱监测。
在对地表环境现象进行连续监测和研究时，常需要获取长时间序列的地表参数。受地表热红外波长较短，无法穿透云雾的影响，当地球表面被云层覆盖时，陆地表面信息存在与云像元混合或完全被云屏蔽的现象。事实上，云-晴空条件占全球每天实际天气的一半以上，这意味着如果在分析过程中排除包含云的像元，我们将失去大量有关陆地表面状况的信息。因此，提出一种对存在云污染的像元的地表温度进行重建的方法，对干旱监测、土壤水分估算及森林火灾等环境现象的连续监测具有重要的意义。
由于LST受太阳辐射、气温、地表覆盖、土壤水分等多种环境因素的影响，LST的重建比反射率的重建更难。因此在云覆盖条件下进行LST重建的研究较少，现有的LST重建算法主要分为两大类，第一类算法主要是重建晴空状态下的LST，但利用邻近时空有效像元重建的LST反映的是该像元在晴空条件下的状态，而非云下地表的真实情况；第二类算法主要是重建云覆盖下地表的真实温度值，这意味着模型中需要充分考虑云层存在的影响，如基于一维热传导方程、日地表温度变化及太阳短波净辐射的计算方法，但如果每天从静止卫星中获得少于6期晴空条件LST数据时该算法将无法使用。
发明内容
本发明要解决的技术问题是提供一种云下像元真实地表温度重建方法，以减少或避免前面所提到的问题。
为解决上述技术问题，本发明提供了一种云下像元真实地表温度重建方法，其通过对连续多日的MODIS日产品数据进行处理，获得指定日期的MODIS日产品数据中的云下像元的真实地表温度，其包括如下步骤：
步骤1，对于待处理的指定日期的MODIS数据，进行像元判断，获得云下像元位置，并将该像元值设定为无效值；同时对与指定日期相邻的多日的MODIS日产品数据进行同样的像元判断。
步骤2，对于步骤1处理后的指定日期的MODIS数据，创建一个9*9像元大小的移动窗口，按照从左至右，从上到下的遍历方式，逐一判断，使窗口中心像元为云下像元，同时创建不同的数组分别存放对应相邻的多日的MODIS日产品数据的移动窗口内的晴空像元值；
步骤3，对于步骤2所创建的移动窗口，确定每个像元对的距离参数(Di)，也即是每个晴空像元相对中心的云下像元的距离参数，Di的计算公式如下：
上式中，(x0，y0)为窗口中心像元在空间上的位置，(xi，yi)为移动窗口内晴空像元所在的位置。
步骤4，计算步骤2的移动窗口中的像元对的环境相似参数Si，计算公式如下：
Si＝Ri×K+1
式中，Si为环境相似参数，K为常数，Ri表征中心像元(x0，y0)和晴空像元(xi，yi)之间的地表环境的差异，其计算公式如下：
Ri＝|Rc-Rv|
Rc＝B(x0，y0，t0)-B(x0，y0，tp)
Rv＝B(xi，yi，t0)-B(xi，yi，tp)
式中，B为地表反射率，Rc为移动窗口内中心像元(x0，y0)地表反射率在邻近时间上的差异；Rv为窗口内晴空像元(xi，yi)地表反射率在邻近时间上的差异。
步骤5，根据步骤3、步骤4中获得的Di和Si，计算每个晴空像元的权重系数Wi，具体计算公式如下：
上式中，N为待处理的指定日期的MODIS数据中的移动窗体内的晴空像元的总个数。
步骤6，计算相关的MODIS日产品数据对于指定日期的MODIS数据的云校正因子CF，具体计算公式如下：
其中，m为与云相关的参数，其计算公式如下：
m＝Rc-Rv
F为用于反映与植被类型、土壤类型及土壤湿度相关的参数，
Gclear为晴空条件下太阳辐射能量，
albedo为地表反照率。
步骤7，重建步骤2的窗口中心的云下像元的真实地表温度，具体按如下公式计算：
上式中，
LST(x0，y0，t0)为移动窗口内中心像元重建后的LST值，
P为数据日期，
N为待处理的指定日期的MODIS数据中的移动窗体内的晴空像元的总个数，
Wi为步骤5中基于距离参数和环境相似参数计算的权重系数，
LST(x0，y0，tp)为中心像元在tp时的LST值，
LST(xi，yi，tp)为邻近晴空像元在tp时的LST值，
LST(xi，yi，t0)为t0时晴空像元LST值，
CFp为步骤6中获得的对应的MODIS的日产品数据计算所得的云校正参数。
优选地，在步骤1中，在计算机上利用IDL读取MOD11A1 LST日产品数据中的QC数据集，将十进制转换为二进制，结合MOD11Collection6用户指导书，判识云下像元，将云下像元的数值设定为无效值0。
优选地，在步骤1中，对与指定日期相邻的八天的MODIS日产品数据进行像元判断。
优选地，在步骤1中，对与指定日期相邻的前四天和后四天的MODIS日产品数据进行像元判断。
优选地，在步骤4中，K的数值设为10，
优选地，在步骤4中，最少只需要获得三组数据。
优选地，在步骤4中，获得八组数据。
优选地，在步骤6中，F设定为140(Wm-2)。
优选地，在步骤7中，Pmin＝t0-4，Pmax＝t0+4。
优选地，在步骤7中，Pmin＝t0-8，Pmax＝t0。
本发明基于云下像元与晴空像元在时间和空间上的联系，提出了一种云下像元真实地表温度重建的方法。根据地理学第一定律，引入距离参数(Di)反映像元对在空间上的相关性；同时结合不同地表类型地表反射率的差异，构建环境相似性参数(Si)用以表征像元对间地表环境的相似性；综合分析像元对在空间分布和地表环境上的差异，通过权重系数(Wi)真实有效地反映晴空像元与云下像元地表温度的相似程度。此外，考虑云层会削弱到达地表的太阳辐射量，从而改变地表温度值，本发明提出云校正参数(CF)表征像元对到达地表辐射量的差异。综合以上参数解决了太阳辐射、地表覆盖、土壤条件以及空间分布等因素对地表温度重建的影响，利用周围晴空像元真实有效地还原了云下像元的地表温度值，为后续多云地区大范围旱情监测、土壤水分估算以及森林火险监测提供了可能。
附图说明
以下附图仅旨在于对本发明做示意性说明和解释，并不限定本发明的范围。其中，
图1为根据本发明的一个具体实施例的一种云下像元真实地表温度重建方法的算法实现的流程原理示意图；
图2为根据图1的流程对云下像元地表温度的重建过程的原理示意图；
图3为图1的实施例所涉及的区域及LST地面站点位置分布图；
图4为图3的LST地面站点所在区域的原始LST值、重建后LST值、地面观测LST值和土壤含水量与时间的关系图；
图5为图4的重建前后LST值与观测值间的误差分析示意图；
图6为图3的重建前后LST遥感图像对比图。
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发明的具体实施方式。
图1为根据本发明的一个具体实施例的一种云下像元真实地表温度重建方法的算法实现的流程原理示意图；图2为根据图1的流程对云下像元地表温度的重建过程的原理示意图；图3为图1的实施例所涉及的区域及LST地面站点位置分布图；图4为图3的LST地面站点所在区域的原始LST值、重建后LST值、地面观测LST值和土壤含水量与时间的关系图；图5为图4的重建前后LST值与观测值间的误差分析示意图；图6为图3的重建前后LST遥感图像对比图。参见图1-6所示，本发明基于云下像元与晴空像元在时间和空间上的联系，提出了一种云下像元真实地表温度重建方法，其包括如下步骤：
步骤1：对于待处理的指定日期的MODIS数据，进行像元判断，获得云下像元位置，并将该像元值设定为无效值0；同时对与指定日期相邻的多日的MODIS日产品数据进行同样的像元判断。
如背景技术所述，MODIS陆地标准数据产品具有不同的时间分辨率，在本发明中，通过处理MODIS的日产品数据来解决云下像元真实地表温度重建的问题，
具体来说，当需要对指定日期的MODIS的日产品数据进行分析处理时，还需要调用与指定日期相邻的多天的MODIS的日产品数据参与计算处理，对于指定日期的MODIS的日产品数据中的云下像元，当该像元在其他时间处于晴空状态时，则可通过本发明所提供的方法来重建获得准确的真实地表温度，理论上来说，相邻的多天的MODIS的日产品数据越多，就越能保障获得指定日期的MODIS的日产品数据中的云下像元在其他日期处于晴空状态。
通常情况下，需要调用八天的MODIS的日产品数据参与对于已发生的特定日期的历史数据的处理，即可满足需要指定日期的MODIS的日产品数据中的云下像元在其他日期有处于晴空的状态的需要，实践中，可调用指定日期前四天和后四天的MODIS的日产品数据，这样最终重建的地表温度准确率即可达到较高的程度，这也是本发明所引用的具体实施例中调用的数据的时间关系。
当然，对于实时的数据的处理，例如，对当日的MODIS的日产品数据进行处理，重建地表温度时，也可调用前八天的MODIS的日产品数据。在本发明涉及的具体实施例中，以指定日期前四天和后四天的MODIS的日产品数据为例进行说明。
在获得待处理的指定日期以及前后四天的MODIS的日产品数据后(也就是共计9天的MODIS的日产品数据)，可在计算机上利用IDL(Interactive Data Language)读取MOD11A1 LST日产品数据中的QC数据集，将十进制转换为二进制，结合MOD11Collection6用户指导书，判识云下像元，也即是云污染及云阴影像元，将云下像元的数值设定为无效值0；
具体来说，例如，对于如图3所示的包含大满站和混合林站的区域，当需要对2018年6月5日的MODIS数据进行处理，重建2018年6月5日的云污染及云阴影像元，也即是云下像元的真实地表温度时，同时也需要调取前四天(也即是2018年6月1日至6月4日)和后四天(也即是2018年6月6日至6月9日)的MODIS数据，在计算机上，对共九天的MODIS数据进行处理，对于每天的MODIS数据，将云下像元的数值设定为无效值0。
步骤2：对于步骤1处理后的数据，创建一个9*9像元大小的移动窗口，按照从左至右，从上到下的遍历方式，逐一判断，使窗口中心像元落在数值为0的像元处(也即是使中心像元为云下像元)，同时创建不同的数组分别存放相邻的多天的移动窗口内的晴空像元值；
如步骤1中所述，在本发明的具体实施例中，调用的是指定日期前四天和后四天的MODIS的日产品数据，因此，在本步骤中，具体来说，为了重建2018年6月5日的云下像元真实地表温度，在计算机上，创建一个9*9像元大小的移动窗口，按照从左至右，从上到下的遍历方式，逐一判断，当窗口中心像元落在数值为0的像元处时，记录移动窗口内晴空像元的位置以及具体的数值，同时需要创建不同的数组分别存放相前后四天(也即是2018年6月1日至6月4日和2018年6月6日至6月9日)的MODIS数据中，对应2018年6月5日当天的晴空像元的位置的像元的具体的数值。
步骤3：利用晴空像元和中心像元(也即是云下像元)的位置坐标，确定每个像元对的距离参数(Di)，也即是每个晴空像元相对中心的云下像元的距离参数，Di的计算公式如下：
上式中，(x0，y0)为窗口中心像元在空间上的位置，(xi，yi)为移动窗口内晴空像元所在的位置。在地表质地均一的条件下，Di值越小，中心像元和晴空像元的地表温度值越接近；反之，像元对之间的地表温度差异越大。
正如背景技术所述，本发明所要处理的是云-晴空条件的MODIS数据，也即是说，在这些数据对应的像元中，除了云下像元，其他都是晴空像元，如步骤2所述，本发明使用的移动窗口大小为9*9像元，这也就意味着(xi，yi)中i的取值范围为-4，-4到4，4。
虽然本发明需要使用共九天的MODIS数据，但因为在这九天中所涉及区域位置是相同的，也就是说对于每天的MODIS数据，在相互关联使用时，9*9像元所覆盖的像元位置相同，因此，在本步骤中，每个晴空像元相对中心的云下像元的距离参数Di只需要计算一次即可。
步骤4：研究表明红光波段反射率与地表土壤水分间存在明显的相关性，土壤水分含量越高，反射率越低；此外，反射率对云的敏感性很高，其常被应用于云检测。而像元间LST的差异常由两像元土壤水分和云量空间差异造成的，两像元间LST差异越小，其地表环境越相似，红光波段发射率的差异也越小。因此，本发明提出利用大气层顶的反射率来计算环境相似参数(Si)，计算公式如下：
Si＝Ri×K+1
式中，Si为环境相似参数，K为常数，用于平衡Si和Di的取值范围，申请人通过大量的数据验证，在本发明中将K的数值设为10，Ri表征中心像元(x0，y0)和晴空像元(xi，yi)之间的地表环境的差异，其计算公式如下：
Ri＝|Rc-Rv|
Rc＝B(x0，y0，t0)-B(x0，y0，tp)
Rv＝B(xi，yi，t0)-B(xi，yi，tp)
式中，B为地表反射率(可从MODIS陆地标准数据产品直接获得)，Rc为移动窗口内中心像元(x0，y0)地表反射率在邻近时间上的差异；Rv为窗口内晴空像元(xi，yi)地表反射率在相邻时间上的差异。由此，Si值越大，表明像元对间地表环境差异越大，像元间地表温度的差异也就越大；反之，像元对间地表温度越接近。
如步骤1中所述，在本发明中，对于指定日期的MODIS的日产品数据中的云下像元，需要该像元在其他时间处于晴空状态，才能获得数据进行后继重建处理，因此，在本发明中，需要调用与指定日期相邻的八天的MODIS的日产品数据，首先分别对每个日产品数据进行中心像元(也即是云下像元)是否处于晴空状态的判断，如中心像元处于晴空状态，则对该对应的日产品数据按照上述公式进行计算，如中心像元同样为云下状态，则该对应的日产品数据不需要按照上述公式进行计算。
也就是说，在调用指定日期前四天和后四天的MODIS的日产品数据的具体实施例中，在本步骤中，Ri和Si共需获得最多八组数据(也就是在所调用的八天的数据中，中心像元均为晴空状态的情况)，发明人通过实践发现，最少只需要有三组数据(也就是在所调用的八天的数据中，中心像元为晴空状态的情况只有三天的数据)，即可重建出满足应用需求(即干旱监测、土壤水分估算及森林火灾等监控需要)的云下像元的真实地表温度，如能有八组数据，则能提高最终重建的真实地表温度的精确度。通常情况下，调用八天的数据即可确保能够获得至少三天的数据。
因此，在上式中，tp用于标识所调用的MODIS的日产品数据与指定日期的MODIS的日产品数据之间的时间关系，具体在本发明的实施例中，当通过前后四天的MODIS的日产品数据进行比对，参与计算时，例如，当(x0，y0，t0)对应2018年6月5日的中心像元时，(x0，y0，tp)的取值则则需要分别对应前四天(也即是2018年6月1日至6月4日)和后四天(也即是2018年6月6日至6月9日)的中心像元数据。也就是说，tp的取值范围为(t0-4，t0+4)。
当然，本领域技术人员应当理解，当需要对实时的数据进行处理，例如，需要对当日的MODIS的日产品数据进行处理，重建地表温度时，也可调用前八天的MODIS的日产品数据，此时，tp的取值范围即为(t0-8，t0)。
步骤5：根据步骤3、步骤4中获得的Di和Si与像元对地表温度间的关系，对于每个晴空像元用权重系数(Wi)来表示该晴空像元在重建中心像元LST值时所占的比重，其计算公式如下：
式中，Di为步骤3中的距离参数(对应每个9*9像元大小的移动窗口，其为唯一数值)，Si为步骤4中的环境相似参数(对应每个9*9像元大小的移动窗口，其可以具有最少三组数值)，N为步骤2中移动窗体内晴空像元的总个数。Wi表示第i个晴空像元在重建中心像元LST值时所占的比重。Wi越大，所占的比重越大，中心像元LST与该像元LST值越相似。
由于在本发明具体实施例中，需要通过前后四天的MODIS的日产品数据进行比对，且在步骤4中，环境相似参数(Si)需要获取最多八组数据，因此，在本步骤中，对于每个9*9像元大小的移动窗口，每个晴空像元同样也会需要获得与步骤4对应的最多八组数据。
步骤6：在太阳辐射能量传输到地表的过程中，云层会削减到达地表的太阳辐射能，从而改变云下像元的LST值。基于能量守恒，本发明提出一个云校正因子用以表征云层削弱的那部分太阳辐射能量：云校正因子CF的具体计算公式如下：
其中，m为与云相关的参数，其计算公式如下：
m＝Rc-Rv
通常云下像元由于有云存在，会使得红光波段反射率要高于晴空像元，因此，对于本发明所处理的MODIS数据来说，由于经由步骤2的定位，已使得9*9像元大小的移动窗口的中心像元落在了云下像元处，如步骤4所述，在本发明中，对于指定日期的MODIS的日产品数据中的云下像元，需要该像元在其他时间处于晴空状态，才能获得数据进行后继重建处理，因此，在本发明中，需要调用与指定日期相邻的八天的MODIS的日产品数据，首先分别对每个日产品数据进行中心像元(也即是云下像元)是否处于晴空状态的判断，如中心像元处于晴空状态，则对该对应的日产品数据按照上述公式进行计算m以及CF，如中心像元同样为云下状态，则该对应的日产品数据不需要按照上述公式进行计算。
也就是说，当(t0、tp)中tp的中心像元为晴空像元时，则其地表反射率的差值(Rc)必然大于t0的晴空像元的地表反射率的差值(Rv)，因此m被定义为Rc-Rv；
F为用于反映与植被类型、土壤类型及土壤湿度相关的参数，申请人通过大量已有数据分析、计算验证，在本发明中将F设定为140(Wm-2)；
Gclear为晴空条件下太阳辐射能量，具体的计算获得方法可参考张春桂等人在《利用卫星资料估算福建晴空太阳辐射》一文中给出的利用MODIS数据计算Gclear方法求得。
albedo为地表反照率，其同样可从MODIS陆地标准数据产品中获得。
步骤7：综合步骤2-6获得的移动窗口内各晴空像元LST值、权重系数和云校正参数，构建云污染像元地表温度重建计算公式，具体按如下公式计算：
上式中，
LST(x0，y0，t0)为移动窗口内中心像元(即云下像元)重建后的LST值，
P为数据日期，在调用指定日期前四天和后四天的MODIS的日产品数据的具体实施例中，Pmin＝t0-4，Pmax＝t0+4，本领域技术人员应当理解，在调用指定日期前八天的MODIS的日产品数据的情况下(即需要对实时的数据进行处理时)，Pmin＝t0-8，Pmax＝t0，
N为步骤5中统计的晴空像元的有效个数，
Wi为步骤5中基于距离参数和环境相似参数计算的权重系数，
LST(x0，y0，tp)为中心像元在tp时的LST值，
LST(xi，yi，tp)为邻近晴空像元在tp时的LST值，
LST(xi，yi，t0)为t0时晴空像元LST值，
CFp为步骤6中获得的对应的MODIS的日产品数据计算所得的云校正参数。
在通过上述步骤1-7完成云下像元真实地表温度重建后，可利用研究区内站点实测的地表温度值与重建后LST间的相关性系数R、偏差Bias以及均方根误差RMSE来评定LST重建的精度。
其中，偏差Bias计算公式如下：
均方根误差RMSE是观测值与真值偏差的平方和观测次数N比值的平方根。其计算公式：
RMSE为均方根误差，N1为观测次数，Xobs，i为重建后的LST值，Xmodel，i为地面站点实测的LST值。
附图3为本实例研究区所在位置分布图，研究区位于中国西北部，地处91.7°E-112.2°E，30.7°N-43.1°N，包括甘肃、宁夏、山西以及蒙古、青海、四川部分地区。图中混合林站和大满站为HiWATER两个地表温度实测站点。
按1～7步骤对研究区内云污染、云阴影像元进行填补重建，图2为云污染像元地表温度值重建模拟过程，基于该算法LST重建后的结果如附图6所示，大量受云层影响的缺失像元LST值得到重建，且在DOY85(儒略日)LST重建结果中能清晰地显示青海湖；附图4中，2014年混合林和大满站点处，实测LST值，原始LST值、重建后LST值以及土壤含水量时间序列上的变化趋势一致，表明基于本发明的重建方法，云下像元LST重建后在时间和空间分布上与地表实际温度保持一致。
表1
通过比较混合林和大满两站点，原始LST和重建后LST与实测的LST间相关性分析结果(表1)，虽然重建后的RMSE相对较大，但重建前后LST与实测的LST间的相关性都较好(R均大于0.90)，且偏差的绝对值介于0～3K之间。附图5为合并混合林和大满站点后，原始LST、重建后LST以及组合原始LST和重建后LST值与实测LST值间的相关性结果，结果显示三组数据与实测LST间的R分别为0.97、0.93和0.96；RMSE分别为3.25K、4.83K和4.21K；偏差分别为1.12K、-1.94K及-0.59K，表明RSDAST算法重建的LST与真实LST间的相关性较高，重建后的LST在一定程度上能较好地表征研究区域某一观测时刻地表的真实温度。
综上所述，基于本发明的重建方法，不仅能实现大面积云下像元LST的重建，且重建后的LST与地面实测的LST间具有较高的相关性，表明该算法适用于多云地区云下像元地表温度的重建，这对多云地区旱情监测和土壤水分估算至关重要。
本领域技术人员应当理解，虽然本发明是按照多个实施例的方式进行描述的，但是并非每个实施例仅包含一个独立的技术方案。说明书中如此叙述仅仅是为了清楚起见，本领域技术人员应当将说明书作为一个整体加以理解，并将各实施例中所涉及的技术方案看作是可以相互组合成不同实施例的方式来理解本实用新型的保护范围。
以上所述仅为本发明示意性的具体实施方式，并非用以限定本发明的范围。任何本领域的技术人员，在不脱离本发明的构思和原则的前提下所作的等同变化、修改与结合，均应属于本发明保护的范围。
Claims (10)
1.一种云下像元真实地表温度重建方法，其特征在于，其通过对连续多日的MODIS日产品数据进行处理，获得指定日期的MODIS日产品数据中的云下像元的真实地表温度，其包括如下步骤：
步骤1，对于待处理的指定日期的MODIS数据，进行像元判断，获得云下像元位置，并将该像元值设定为无效值；同时对与指定日期相邻的多日的MODIS日产品数据进行同样的像元判断，
步骤2，对于步骤1处理后的指定日期的MODIS数据，创建一个9*9像元大小的移动窗口，按照从左至右，从上到下的遍历方式，逐一判断，使窗口中心像元为云下像元，同时创建不同的数组分别存放对应相邻的多日的MODIS日产品数据的移动窗口内的晴空像元值；
步骤3，对于步骤2所创建的移动窗口，确定每个像元对的距离参数（
），也即是每个晴空像元相对中心的云下像元的距离参数，
的计算公式如下：
上式中，
为窗口中心像元在空间上的位置，（
）为移动窗口内晴空像元所在的位置，
步骤4，计算步骤2的移动窗口中的像元对的环境相似参数
，计算公式如下：
式中，
为环境相似参数，K为常数，
表征中心像元
和晴空像元（
）之间的地表环境的差异，其计算公式如下：
式中，B为地表反射率，
为移动窗口内中心像元
地表反射率在邻近时间上的差异；
为窗口内晴空像元（
）地表反射率在邻近时间上的差异，
步骤5，根据步骤3、步骤4中获得的
和
，计算每个晴空像元的权重系数
，具体计算公式如下：
上式中，N为待处理的指定日期的MODIS数据中的移动窗体内的晴空像元的总个数，
步骤6，计算相关的MODIS日产品数据对于指定日期的MODIS数据的云校正因子CF，具体计算公式如下：
其中，
为与云相关的参数，其计算公式如下：
F为用于反映与植被类型、土壤类型及土壤湿度相关的参数，
步骤7，重建步骤2的窗口中心的云下像元的真实地表温度，具体按如下公式计算：
上式中，
P为数据日期，
N为待处理的指定日期的MODIS数据中的移动窗体内的晴空像元的总个数，
2.根据权利要求1所述的方法，其特征在于，在步骤1中，在计算机上利用IDL读取MOD11A1 LST日产品数据中的QC数据集，将十进制转换为二进制，结合MOD11Collection6用户指导书，判识云下像元，将云下像元的数值设定为无效值0。
3.根据权利要求1所述的方法，其特征在于，在步骤1中，对与指定日期相邻的八天的MODIS日产品数据进行像元判断。
4.根据权利要求3所述的方法，其特征在于，在步骤1中，对与指定日期相邻的前四天和后四天的MODIS日产品数据进行像元判断。
5.根据权利要求1所述的方法，其特征在于，在步骤4中，K的数值设为10。
6.根据权利要求1所述的方法，其特征在于，在步骤4中，最少只需要获得三组数据。
7.根据权利要求3所述的方法，其特征在于，在步骤4中，获得八组数据。
8.根据权利要求1所述的方法，其特征在于，在步骤6中，F设定为140(
)。
9.根据权利要求1所述的方法，其特征在于，在步骤7中，Pmin=
，Pmax=
。
10.根据权利要求1所述的方法，其特征在于，在步骤7中，Pmin=
，Pmax=
。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910609171.7A CN110348107B (zh) | 2019-07-08 | 2019-07-08 | 一种云下像元真实地表温度重建方法 |
| AU2020101048A AU2020101048A4 (en) | 2019-07-08 | 2020-06-18 | Method for reconstructing true land surface temperature of cloudy pixel |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910609171.7A CN110348107B (zh) | 2019-07-08 | 2019-07-08 | 一种云下像元真实地表温度重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110348107A CN110348107A (zh) | 2019-10-18 |
| CN110348107B true CN110348107B (zh) | 2022-06-24 |
Family
ID=68178215
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910609171.7A Expired - Fee Related CN110348107B (zh) | 2019-07-08 | 2019-07-08 | 一种云下像元真实地表温度重建方法 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN110348107B (zh) |
| AU (1) | AU2020101048A4 (zh) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111915694B (zh) * | 2020-07-01 | 2022-11-08 | 河海大学 | 一种顾及时空特征的云覆盖像元地表温度重建方法 |
| CN114019579B (zh) * | 2021-08-24 | 2023-11-24 | 中国农业科学院农业资源与农业区划研究所 | 高时空分辨率近地表空气温度重构方法、系统、设备 |
| CN114034651A (zh) * | 2021-11-10 | 2022-02-11 | 北京环境特性研究所 | 一种全球地表光谱基础数据的生成方法和装置 |
| CN114218756A (zh) * | 2021-11-24 | 2022-03-22 | 中国农业科学院农业资源与农业区划研究所 | 一种基于地表温度年变化模型的云下地表温度重建方法 |
| CN114913296B (zh) * | 2022-05-07 | 2023-08-11 | 中国石油大学(华东) | 一种modis地表温度数据产品重建方法 |
| CN115014579A (zh) * | 2022-06-01 | 2022-09-06 | 联通(四川)产业互联网有限公司 | 一种基于大数据的农业环境温度监测方法及系统 |
| CN116011342B (zh) * | 2023-02-08 | 2023-07-21 | 中国自然资源航空物探遥感中心 | 一种高分辨率热红外地表温度全天候重建方法 |
| CN118411622B (zh) * | 2024-04-30 | 2024-10-18 | 西南交通大学 | 一种30米分辨率全天空地表温度估算方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105841819A (zh) * | 2016-03-23 | 2016-08-10 | 中国科学院遥感与数字地球研究所 | 一种有云条件下的地表温度的估算方法和装置 |
| CN106169058A (zh) * | 2015-12-25 | 2016-11-30 | 中国石油大学(华东) | 一种基于微波遥感与时空信息的云下像元lst估算方法 |
| WO2018145229A1 (zh) * | 2017-02-10 | 2018-08-16 | 广西壮族自治区气象减灾研究所 | 一种近地面气温的大面积精确反演方法 |
| CN108984867A (zh) * | 2018-06-28 | 2018-12-11 | 中国科学院地理科学与资源研究所 | 基于gf-4和modis结合的蒸散发遥感反演方法和系统 |
| CN109186774A (zh) * | 2018-08-30 | 2019-01-11 | 清华大学 | 地表温度信息获取方法、装置、计算机设备和存储介质 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4365887B1 (ja) * | 2009-03-18 | 2009-11-18 | 株式会社パスコ | 地表面画像データの生成方法および生成装置 |
-
2019
- 2019-07-08 CN CN201910609171.7A patent/CN110348107B/zh not_active Expired - Fee Related
-
2020
- 2020-06-18 AU AU2020101048A patent/AU2020101048A4/en not_active Ceased
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106169058A (zh) * | 2015-12-25 | 2016-11-30 | 中国石油大学(华东) | 一种基于微波遥感与时空信息的云下像元lst估算方法 |
| CN105841819A (zh) * | 2016-03-23 | 2016-08-10 | 中国科学院遥感与数字地球研究所 | 一种有云条件下的地表温度的估算方法和装置 |
| WO2018145229A1 (zh) * | 2017-02-10 | 2018-08-16 | 广西壮族自治区气象减灾研究所 | 一种近地面气温的大面积精确反演方法 |
| CN108984867A (zh) * | 2018-06-28 | 2018-12-11 | 中国科学院地理科学与资源研究所 | 基于gf-4和modis结合的蒸散发遥感反演方法和系统 |
| CN109186774A (zh) * | 2018-08-30 | 2019-01-11 | 清华大学 | 地表温度信息获取方法、装置、计算机设备和存储介质 |
Non-Patent Citations (2)
| Title |
|---|
| Reconstructing daily clear-sky land surface temperature for cloudy regions from MODIS data;Liang Sun 等;《Computers and Geosciences》;20170831(第105期);全文 * |
| 江苏地区MODIS LST产品重建研究;严婧等;《环境科学与技术》;20140115(第01期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110348107A (zh) | 2019-10-18 |
| AU2020101048A4 (en) | 2020-07-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110348107B (zh) | 一种云下像元真实地表温度重建方法 | |
| Wang et al. | Intercomparisons of cloud mask products among Fengyun-4A, Himawari-8, and MODIS | |
| Zhang et al. | An operational approach for generating the global land surface downward shortwave radiation product from MODIS data | |
| Wan et al. | Validation of the land-surface temperature products retrieved from Terra Moderate Resolution Imaging Spectroradiometer data | |
| Ye et al. | Land surface temperature estimate from Chinese Gaofen-5 satellite data using split-window algorithm | |
| Ermida et al. | Validation of remotely sensed surface temperature over an oak woodland landscape—The problem of viewing and illumination geometries | |
| Townshend et al. | The 1 km resolution global data set: needs of the International Geosphere Biosphere Programme | |
| Leroy et al. | Sun and view angle corrections on reflectances derived from NOAA/AVHRR data | |
| Vermote et al. | Atmospheric correction algorithm: spectral reflectances (MOD09) | |
| Jia et al. | Generating a 2-km, all-sky, hourly land surface temperature product from Advanced Baseline Imager data | |
| Chen et al. | A data-model fusion approach for upscaling gross ecosystem productivity to the landscape scale based on remote sensing and flux footprint modelling | |
| Li et al. | Atmospheric correction of geostationary satellite ocean color data under high solar zenith angles in open oceans | |
| CN102298150B (zh) | 全球陆表宽波段发射率反演方法及系统 | |
| Carrer et al. | Land surface albedo derived on a ten daily basis from Meteosat Second Generation Observations: The NRT and climate data record collections from the EUMETSAT LSA SAF | |
| Sekertekin et al. | Modeling diurnal Land Surface Temperature on a local scale of an arid environment using artificial Neural Network (ANN) and time series of Landsat-8 derived spectral indexes | |
| Carrer et al. | Comparing operational MSG/SEVIRI land surface albedo products from Land SAF with ground measurements and MODIS | |
| Leng et al. | A method for deriving all‐sky evapotranspiration from the synergistic use of remotely sensed images and meteorological data | |
| Varma | Measurement of precipitation from satellite radiometers (visible, infrared, and microwave): Physical basis, methods, and limitations | |
| Yu et al. | All-sky total and direct surface shortwave downward radiation (SWDR) estimation from satellite: Applications to MODIS and Himawari-8 | |
| Wan et al. | Accuracy Evaluation and Parameter Analysis of Land Surface Temperature Inversion Algorithm for Landsat‐8 Data | |
| Xu et al. | Reconstructing all-weather daytime land surface temperature based on energy balance considering the cloud radiative effect | |
| Wu et al. | Improving the capability of water vapor retrieval from Landsat 8 using ensemble machine learning | |
| Zhao et al. | A Nonlinear Split-Window Algorithm for Retrieving Land Surface Temperatures from Fengyun-4B Thermal Infrared Data | |
| Ma et al. | An in-flight radiometric calibration method considering adjacency effects for high-resolution optical sensors over artificial targets | |
| Huang et al. | A 3D approach to reconstruct continuous optical images using lidar and MODIS |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220624 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |