一种获取全天候有效土壤水分的方法

摘要

本发明公开了一种获取全天候有效土壤水分的方法，包括以下步骤：采集CLDAS气象数据；采集MODIS数据，利用ENVI软件进行预处理，最后输出为tiff格式；对MODIS像云进行判别，分为MODIS无云像元与MODIS有云像元，基于CLDAS气象数据和MODIS数据进行估算MODIS无云像元和MODIS有云像元的有效土壤水分。本发在现有遥感反演有效土壤水分的基础上，针对有云像元无法正常反演有效土壤水分的情况，提出利用气象数据信息进行估算的方法，从而实现获取全天候的有效土壤水分数据。

权利要求

1.一种获取全天候有效土壤水分的方法，其特征在于，包括以下步骤：步骤1、采集CLDAS气象数据：CLDAS气象数据包括气温、风速、短波辐射和比湿；

步骤2、采集MODIS数据，利用ENVI软件进行预处理，最后输出为tiff格式，作为输入数据：MODIS数据包括地表温度产品、8天合成的叶面积指数产品以及16天合成归一化植被指数产品；

步骤3、对MODIS数据的像元进行判别，分为MODIS无云像元与MODIS有云像元，基于CLDAS气象数据和MODIS数据进行估算MODIS无云像元和MODIS有云像元的有效土壤水分；

其中：步骤3中的估算MODIS无云像元的有效土壤水分具体为：计算干燥裸土温度Ts,max，湿润裸土温度Ts,min，受水分胁迫的全植被覆盖温度Tc,max和水分状况良好的全植被覆盖温度Tc,min四个温度值，计算公式分别为：其中 ，Δ是 饱 和水 汽 压随 温 度变 化的 斜 率 ，表 示 为气 温的 函 数 ；

其中，Ta是气温，从CLDAS数据

中读取；

Cv是空气比热，取值为1295.16；

rcx和rcp分别是最大冠层阻抗和潜在蒸发时的冠层阻抗，rcx取值为2000，rcp根据LAI进行计算，rcp＝100/LAI，其中，LAI是叶面积指数，LAI从MODIS数据中读取；

γ是干湿表常数，取值0.066；

VPD是饱和水汽压差，利用CLDAS的气温和比湿数据进行计算；

Rn是净辐射，利用CLDAS数据读取下行短波辐射Sd、气温Ta以及四个极端条件下的反照率进行计算；Rn＝(1-αs)Sd+εsεaσTa4-εsσTa4其中，αs是地表反照率，在干燥裸土、湿润裸土、受水分胁迫的全植被覆盖和水分良好的全植被覆盖条件下分别取值为0.20，0.10，0.25和

0.15；εs为地表比辐射率，干燥和湿润的裸土条件分别为0.96和0.97，受水分胁迫的全植被和水分良好的全植被分别为0.975和0.985；σ是Stefan-Boltzmann常数，σ＝5.67×10-8；εa是天空发射率，表示为气温的函数εa＝9.2×10-6×Ta2；

G是土壤热通量，表示为Rn的函数；G＝[0.05+0.265×(1-FVC)]Rn，其中，FVC是植被覆盖度，利用归一化植被指数NDVI产品MOD13A2计算：其中，NDVI直接从MODIS数据中读取，NDVImax和NDVImin分别对应全植被覆盖和裸土的NDVI，分别给定为0.86和0.2；

ra是空气动力学粗糙度，计算公式为:

其中，Z0M和Z0H分别是动量传输和能量传输粗糙度长度，d是零平面位移，z0是表面粗糙长度，u2是2m处风速；设植被高度为h，d＝0.65h，z0＝0.1h，Z0M＝0.1h，Z0H＝0.1Z0M；2m处风速u2由CLDAS数据的10m处风速U进行换算，计算公式为：对于MODIS无云像元，有效土壤水分的计算公式为：其中，Ts是直接从MODIS的温度产品中读取的数据；

再者：步骤3中的估算MODIS有云像元的有效土壤水分具体为：计算地表阻抗：

其中，us是地表粗糙度影响最小的高度处的风速，利用叶面积指数，叶子直径和植被高度进行估算，设植被高度为h，叶子直径为s，us为：其中，LAI是叶面积指数，从MODIS数据直接读取；

MODIS有云像元有效土壤水分M0近似计算为：

说明书

一种获取全天候有效土壤水分的方法

技术领域

[0001] 本发明属于定量遥感技术领域，具体地说，涉及一种获取全天候有效土壤水分的方法。

背景技术

[0002] 有效土壤水分是土壤中能够被作物所利用的水分，它对于作物生长非常重要，也是进行干旱监测的重要指标。在区域尺度中，通常利用光学与热红外遥感技术来获取区域连续的有效土壤水分。然而，当遇到雨云等不利天气条件时，光学与热红外像元可能无法获取有效的观测值，导致遥感方法反演失效，从而无法获得全天候的有效土壤水分。对此，本发明充分挖掘气象数据的潜力，利用现有全天候气象数据产品信息弥补遥感数据缺失给全天候有效土壤水分反演带来的不便，提出一种获取全天候有效土壤水分的方法，为指导农业生产活动和作物生长提供可靠的土壤水分信息，就有重要的应用价值和实际意义。

[0003] 传统的区域有效土壤水分一般利用光学与热红外遥感技术获取，只能在晴空条件下使用，对于非晴空像元，一般无法获取有效的观测值，从而无法反演全天候有效土壤水分。

发明内容

[0004] 有鉴于此，本发明针对现有光学与热红外遥感技术中存在的非晴空像元不能够直接获取有效土壤水分的问题，提供了一种获取全天候有效土壤水分的方法，从而实现全天候有效土壤水分的获取。

[0005] 为了解决上述技术问题，本发明公开了一种获取全天候有效土壤水分的方法，包括以下步骤：

[0006] 步骤1、采集CLDAS气象数据：CLDAS气象数据包括气温、风速、短波辐射和比湿；

[0007] 步骤2、采集MODIS数据，利用ENVI软件进行预处理，最后输出为tiff格式，作为输入数据：MODIS数据包括地表温度产品、8天合成的叶面积指数产品以及16天合成归一化植被指数产品；

[0008] 步骤3、对MODIS数据进行判别，分为MODIS无云像元与MODIS有云像元，基于CLDAS气象数据和MODIS数据进行估算MODIS无云像元和MODIS有云像元的有效土壤水分。

[0009] 进一步地，MODIS无云像元与MODIS有云像元是MODIS产品数据自带的属性，只要下载数据后就能直观获得此信息。

[0010] 进一步地，步骤3中的估算MODIS无云像元的有效土壤水分具体为：

[0011] 计算干燥裸土温度Ts,max，湿润裸土温度Ts,min，受水分胁迫的全植被覆盖温度Tc,max和水分状况良好的全植被覆盖温度Tc,min四个温度值，计算公式分别为：

[0012]

[0013]

[0014]

[0015]

[0016] 其中，Δ是饱和水汽压随温度变化的斜率，表示为气温的函数；其中，Ta是气温，从CLDAS数据
中读取；

[0017] Cv是空气比热，取值为1295.16；

[0018] rcx和rcp分别是最大冠层阻抗和潜在蒸发时的冠层阻抗，rcx取值为2000，rcp根据LAI进行计算，rcp＝100/LAI，其中，LAI从MODIS数据中读取；

[0019] γ是干湿表常数，取值0.066；

[0020] VPD是饱和水汽压差，利用CLDAS的气温和比湿数据进行计算；

[0021] Rn是净辐射，利用CLDAS数据读取下行短波辐射(Sd)、气温(Ta)以及四个极端条件下的反照率进行计算；Rn＝(1-αs)Sd+εsεaσTa4-εsσTa4其中，αs是地表反照率，在干燥裸土、湿润裸土、受水分胁迫的全植被覆盖和水分良好的全植被覆盖条件下分别取值为0.20，0.10，0.25和0.15；εs为地表比辐射率，干燥和湿润的裸土条件分别为0.96和0.97，受水分胁迫的全植被和水分良好的全植被分别为0.975和0.985；σ是Stefan-Boltzmann常数，σ＝5.67×
10-8；εa是天空发射率，表示为气温的函数εa＝9.2×10-6×Ta2；

[0022] G是土壤热通量，表示为Rn的函数；G＝[0.05+0.265×(1-FVC)]Rn，其中，[0023] FVC是植被覆盖度，利用归一化植被指数NDVI产品MOD13A2计算：

[0024]

[0025] 其中，NDVI直接从MODIS数据中读取，NDVImax和NDVImin分别对应全植被覆盖和裸土的NDVI，从MOD13A2数据中得到，或者根据经验公式得到，分别给定为0.86和0.2；

[0026] ra是空气动力学粗糙度，计算公式为:

[0027]

[0028] 其中，Z0M和Z0H分别是动量传输和能量传输粗糙度长度，d是零平面位移，z0是表面粗糙长度，u2是2m处风速；设植被高度为h，d＝0.65h，z0＝0.1h，Z0M＝0.1h，Z0H＝0.1Z0M；2m处风速u2由CLDAS数据的10m处风速U进行换算，计算公式为：

[0029]

[0030] 对于MODIS无云像元，有效土壤水分的计算公式为：

[0031]

[0032] 其中，Ts是直接从MODIS的温度产品中读取的数据。

[0033] 进一步地，步骤3中的估算MODIS有云像元的有效土壤水分具体为：

[0034] 计算地表阻抗：

[0035]

[0036] 其中，us是地表粗糙度影响最小的高度处的风速，利用叶面积指数，叶子直径和植被高度进行估算，设植被高度为h，叶子直径为s，us为：

[0037]

[0038] 其中，LAI是叶面积指数，从MODIS数据直接读取；

[0039] MODIS有云像元有效土壤水分M0近似计算为：

[0040]

[0041] 与现有技术相比，本发明可以获得包括以下技术效果：

[0042] 1)本发明在现有遥感反演有效土壤水分的基础上，针对有云像元无法正常反演有效土壤水分的情况，提出利用气象数据信息进行估算的方法，弥补了单纯使用遥感方法在有云条件下无法反演有效土壤水分的缺陷。

[0043] 2)本发明融合遥感和气象数据信息，能够获取空间上连续的全天候有效土壤水分，为区域研究提供完整的土壤水分数据；

[0044] 3)本发明满足了定量遥感领域研究对地面数据信息完整性的需求。

[0045] 当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有技术效果。

附图说明

[0046] 此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

[0047] 图1是本发明获取全天候有效土壤水分的方法的流程图；

[0048] 图2是本发明2015年9月20日河南省MODIS地表温度；

[0049] 图3是本发明2015年9月20日河南省MODIS归一化植被指数(NDVI)；

[0050] 图4是本发明2015年9月20日河南省气象数据，其中，(a)为比湿，(b)为风速；

[0051] 图5是2015年9月20日河南省有效土壤水分，其中，(a)为传统遥感反演方法，(b)为本发明的方法。

具体实施方式

[0052] 以下将配合实施例来详细说明本发明的实施方式，藉此对本发明如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。

[0053] 本发明获取全天候有效土壤水分所需的数据为中国气象局陆面数据同化系统(China Meteorological Administration Land Data Assimilation System，CLDAS)实时气象数据和美国的MODIS(Moderate Resolution Imaging Spectroradiometer)数据，无需其它辅助数据。

[0054] 本发明提供一种获取全天候有效土壤水分的方法，如图1所示，包括以下步骤：

[0055] 步骤1：CLDAS气象数据：

[0056] 从中国气象数据共享网下载CLDAS气象要素产品，包括气温、风速、短波辐射和比湿，下载网址为：http://data.cma.cn/data/detail/dataCode/NAFP_CLDAS_RT.html。无需进行任何预处理操作，可直接作为输入数据。

[0057] 步骤2：MODIS数据

[0058] 从MODIS数据网站https://ladsweb.nascom.nasa.gov/search/下载产品数据，包括地表温度LST产品MOD11A1、8天合成的叶面积指数(LAI)产品(M OD15A2)以及16天合成归一化植被指数NDVI产品MOD13A2。利用ENVI软件进行预处理，主要包括：投影、拼接、裁剪等，最后输出为tiff格式，作为本发明的输入数据。

[0059] 步骤3：估算MODIS像元有效土壤水分

[0060] 对MODIS像云进行判别，分为MODIS无云像元与MODIS有云像元，MODIS无云像云一般有正常的像元温度值，可以直接利用遥感数据反演有效土壤水分；MODIS有云像元则没有像元温度值，因而不能直接反演有效土壤水分。这是MODIS产品数据自带的属性，只要下载数据后就能直观获得此信息，无需额外方法和过程进行判断；本发明根据不同的情况进行处理：

[0061] (1)无云象元

[0062] 对于MODIS无云像元，首先计算四个温度值，即：干燥裸土温度Ts,max，湿润裸土温度Ts,min，受水分胁迫的全植被覆盖温度Tc,max和水分状况良好的全植被覆盖温度Tc,min，计算公式分别为：

[0063]

[0064]

[0065]

[0066]

[0067] 其中，Δ是饱和水汽压随温度变化的斜率，可以表示为气温的函数，

[0068] 其中，Ta是气温；Cv是空气比热(取值1295.16)；rcx和rcp分别是最大冠层阻抗和潜在蒸发时的冠层阻抗，rcx可以给一定值(2000)，rcp可以根据LAI进行计算(rcp＝100/LAI)，其中，LAI可以从MODIS数据中读取；γ是干湿表常数，取值0.066；Ta是气温，从CLDAS数据中读取；VPD是饱和水汽压差，可以利用CLDAS的气温和比湿数据进行计算；Rn是净辐射，利用CLDAS数据读取下行短波辐射(Sd)、气温(Ta)以及四个极端条件下的反照率进行计算，Rn＝(1-αs)Sd+εsεaσTa4-εsσTa4，其中，αs是地表反照率，在干燥裸土、湿润裸土、受水分胁迫的全植被覆盖和水分良好的全植被覆盖条件下分别取值为0.20，0.10，0.25和0.15。εs为地表比辐射率，干燥和湿润的裸土条件分别为0.96和0.97，受水分胁迫的全植被和水分良好的全植被分别为0.975和0.985。σ是Stefan-Boltzmann常数，σ＝5.67×10-8。εa是天空发射率，表-6 2示为气温的函数εa＝9.2×10 ×Ta 。G是土壤热通量，可以表示为Rn的函数。G＝[0.05+
0.265×(1-FVC)]Rn，其中，FVC是植被覆盖度，可以利用归一化植被指数(NDVI)产品(MOD13A2)计算：

[0069] 其中，NDVI可直接从MODIS数据中读取，NDVImax和NDVImin分别对应全植被覆盖和裸土的NDVI，可以从MOD13A2数据中得到，也可以根据经验公式，分别给定为0.86和0.2。

[0070] ra是空气动力学粗糙度，计算公式为:

[0071]

[0072] 其中，Z0M和Z0H分别是动量传输和能量传输粗糙度长度，d是零平面位移，z0是表面粗糙长度，u2是2m处风速。假设植被高度为h，d＝0.65h，z0＝0.1h，Z0M＝0.1h，Z0H＝0.1Z0M；2m处风速u2由CLDAS数据的10m处风速U进行换算，计算公式为：

[0073]

[0074] 对于MODIS无云像元，有效土壤水分的计算公式为：

[0075]

[0076] 其中，Ts是直接从MODIS的温度产品中读取的数据。

[0077] (2)有云象元

[0078] 对于有云像元，首先计算地表阻抗：

[0079]

[0080] 其中，us是地表粗糙度影响最小的高度处的风速，可以利用叶面积指数，叶子直径和植被高度进行估算。植被高度为h，叶子直径为s，us为：

[0081]

[0082] 其中，LAI是叶面积指数，可以从MODIS数据直接读取。

[0083] MODIS有云像元有效土壤水分M0可以近似计算为：

[0084]

[0085] 根据公式(8)和公式(11)，分别可以计算有云像元和无云像元的有效土壤水分，从而实现全天候有效土壤水分的获取。

[0086] 实施例1

[0087] 以河南省为研究区，图2是2015年9月20日河南省MODIS的地表温度，从图中可以看出河南北部、东部部分地区以及西部部分地区数据缺失，呈现大量空白区域；图3为归一化植被指数。这两个数据是本发明的直接输入数据。图4是本发明的另外一套直接输入数据，即CLDAS气象数据。

[0088] 图5是有效土壤水分反演结果。图5(a)是传统遥感方法反演的结果，图5(b)是本发明的方法。从两个图的对比可以看出，传统遥感反演方法在无有效MODIS像元值的情况下无法反演有效土壤水分，而本发明的方法能够弥补这一缺憾，实现有效土壤水分的全天候反演，且反演结果呈现良好的空间连续性。

[0089] 上述说明示出并描述了发明的若干优选实施例，但如前所述，应当理解发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。