一种确定SM监测中LST/FVC空间的方法

摘要

本申请提供了一种确定SM监测中LST/FVC空间的方法，该方法引入了Biome‑BGC模型中平行阻抗的方式，计算干边的温度端元，简化了干边温度的计算。此外，该方法还发展了一种确定湿边温度端元的逐步选择法。通过计算土壤湿度指数，比较了本方法与传统的理论计算方法。SAMPVEX12综合实验数据、MODIS遥感数据、以及NLDAS‑2气象数据被用于实验研究，结果表明在监测土壤湿度的空间变化上，本申请的方法与之前存在的理论计算方法具有非常高的一致性。

说明书

技术领域

本发明涉及一种确定SM监测中LST/FVC空间的方法。

背景技术

土壤湿度(SM)是陆地表面蒸散(ET)和植被/作物生长的重要水源。它通过控制可利用能量对潜热和显热通量的分配，以及降水对地表径流和入渗的分配，在地表水和能量循环中发挥着重要作用。因此，它是广泛应用的关键参数，如农业生产(Dobriyal等人2012)，干旱监测和预测(AghaKouchak等人2015；Sun等人2013)，水资源管理(Robinson等人，2008)，天气预报(Dai，Trenberth和Qian2004)，以及气候变化(Anderson等人2007；Petropoulos，Ireland和Barrett 2015)等。

遥感是一种在一定时间间隔内从区域尺度到全球尺度获取空间分布的技术。基于遥感数据检索SM的方法几乎涵盖了从可见波段到微波波段的整个光谱范围(Petropoulos、Ireland和Barrett，2015年)。其中，基于微波波段特别是L波段的方法被认为是获得SM的最佳选择，因为波长更长并且与SM更直接相关(Petropoulos，Ireland和Barrett 2015；Peng等2017；Colliander等2017；Wu等2014；Administration 2014)。然而，基于可见光、近红外和热红外波段的方法仍然很重要，因为它们通常具有较长的历史，具有高时空分辨率和丰富的数据源。此外，光学/热遥感数据是将微波SM数据空间分辨率提高到几公里甚至几十米的重要来源(Peng等人2017)。

遥感地表温度和植被盖度或光谱植被指数(LST/FVC)特征空间是用于从光学/热遥感监测SM的概念模型。它可以简单地表示为二维散点图，其中水平轴是FVC从0到1变化，垂直轴是从潜在最小值到潜在最大值的LST(Land Surface Temperature，地表温度)(Carlson 2007；Petropoulos等人2009；Sun2016)。潜在的最小值和最大值对应于陆地表面ET的两个极端条件，其中前者处于饱和水供应而后者处于最大水压力下。所有具有最大水应力的像素形成LST/FVC空间的上边界(即干边)，而具有饱和水供应的像素形成LST/FVC空间的下边界(即湿边)(Sun等人2017)。像素相对于LST/FVC空间中干边和湿边的位置是监测SM和进一步ET的变化的有力指标。例如，温度植被干旱指数(TVDI)(Sandholt、Rasmussen和Andersen 2002)和植被温度条件指数(VTCI)(Wang等人2001)是基于LST/FVC空间监测SM的两个典型指标。此外，它对微波遥感观测估算的降尺度粗分辨率SM也起着重要作用。例如，Merlin等人，提出了一种物理降尺度方法，即基于物理和理论尺度的分解(DISPATCH)，其中LST/FVC空间用于计算显着的降尺度因子土壤蒸发效率(SEE)(Merlin等2012；Merlin等2013)。

将LST/FVC空间应用于SM监测的一个关键问题是确定干边和湿边。视觉识别和自动拟合是直接从遥感图像本身确定干湿边的两种简单方法。然而，它们对研究区域的要求非常苛刻，例如在该区域必须存在干湿边条件下的像素，这通常是无法满足(Tang、Li和Tang 2010)。而且它是经验方法，具有经验方法的所有特征。例如；依赖于场景，缺乏对其他区域的可转移性，不能微调以及描述物理过程。除上述经验方法外，Moran等(1994)，Long等(2012)，Sun等(2016年、2017年)还有三种典型的理论计算方法。

然而，理论计算方法需要相对复杂的迭代过程和参数，例如确定空气动力学阻抗需要的风速和植被高度，这些参数不容易获取(Long，Singh和Scanlon2012；Sun 2016b；Sun等2017；Moran等1994年；Sun等2012)。

发明内容

鉴于上述问题，本申请目的是提供一种确定SM监测中LST/FVC空间的方法，以供SM监测的操作性和非迭代性使用。

本发明的一种确定SM监测中LST/FVC空间的方法，其通过确定LST/FVC空间的湿边的温度端元和干边的温度端元来确定LST/FVC空间；

其中，湿边的温度端元通过式2确定：

其中LST_w代表土地覆盖类型“永久湿地”或“水体”的地表温度，单位K；LST_{no_ice}代表没有“冰雪”的土地覆盖类型以及上述湿地和水体上的地表温度，单位K；T_a是实时近地面气温，单位K；当LST_w存在于一个局部研究区域时，LST_w的平均值被认为是T_s^min和T_c^min；如果不存在LST_w且存在LST_{no_ice}<T_a，则将LST_{no_ice}的最小值设置为T_s^min和T_c^min；对于其他情况，T_a被设为T_s^min。

优选地，干边上的温度端元T_s^max和T_c^max可以通过式3、4、5确定：

其中α_c和α_s是植被组分和裸土组分短波的反照率；ε_c和ε_s是植被组分和裸土组分的宽波段发射率；r_a^c是干燥植被冠层的空气动力学阻抗，是干燥裸土组分的空气动力学阻抗；n_s是分数系数，为0.35；S_d是下行短波辐射；σ是Stefan-Boltzmann常数；ρ是空气密度；c_p是空气比热；ε_a是大气发射率；

rh是植被表面的对流传热阻抗；rr是辐射传热阻抗；gl_sh为叶尺度边界层导度，其值根据土地覆盖类型确定；

r_hs是裸土表面的对流传热阻抗；rbl_min表示蒸汽输送的总空气动力阻力，其值设定为95.0s/m。

优选地，ε_c为0.983，ε_s的值为0.959。

优选地，ε_a＝1.0-0.261exp[-7.77×10^-4(273-T_a)²]。

本申请的确定SM监测中LST/FVC空间的方法，引入了Biome-BGC模型中平行阻抗的方式，计算干边的温度端元，简化了干边温度的计算。此外，该方法还发展了一种确定湿边温度端元的逐步选择法。

附图说明

图1为LST/FVC空间中干边和湿边的示意图；

图2为6月27日的土壤湿度和各种方法计算的SMI；

图3为7月5日的土壤湿度和各种SMI；

图4为7月13日的土壤湿度和各种SMI；

图5为7月17日的土壤湿度和各种SMI；

图6为7月13日的气温、地表温度、地表温度与气温差异以及差异分布；

图7为月27日，7月5日，7月13日，7月17日四天的土壤湿度和各种SMI比较的散点图；

图8为所有时间点上的相关系数随土壤湿度平均值的变化情况；

图9为温度极值变化与输入变量的相对变化。

具体实施方式

下面结合附图对本申请进行详细说明。

LST/FVC空间干湿边界的定义

图1显示了LST/FVC空间中干边和湿边的示意图。目前，对LST/FVC空间的理论形式有不同的解释，即梯形(Long和Singh 2012；Yang和Shang 2013)和两级段梯形(Sun 2016)。但这些方法干边和湿边的解释是相似的，点A和D分别代表具有最大水分胁迫的土壤和植被组分T_s^max和T_c^max，B点和C点代表饱和水供应的土壤和植被组分T_s^min和T_c^min。上述点A、B、C和D也被称为LST/FVC空间的温度端元，线BC和线AD即为湿边和干边。湿边BC代表ET具有饱和供水的地表状态，相应地，干边AD代表ET具有最大水分胁迫的地表状态。

半经验算法确定干湿边

参考之前的研究，干边和湿边通常被认为是线性函数(Petropoulos等2009；卡尔森2007年)，可表示为：

式中，LST_dry和LST_wet代表干边和湿边上的地表温度，随f_v(即FVC)变化，T_s^max,T_s^min,T_c^min,T_c^max对应前面提到的四个端点温度分别对应于极干燥的土壤、极潮湿的土壤、极潮湿的植被和极干燥的植被。

本申请的方法是基于半经验算法而提出的。

本申请中设计了一种逐步选择的方式确定湿边上的温度端元。具体地，如果在一个局地范围内存在水体和永久性湿地，那么该区域上的地表温度平均值，被确定为湿边上的温度端元。如果不存在水体或永久性湿地，那么比较其他地表覆盖类型(不包括冰雪覆盖)上的地表温度与近地表气温。当存在地表温度小于近地表气温的情况时，该局地范围内地表温度的最小值确定为湿边温度端元。当不存地表温度小于近地表气温的情况时，取近地表气温作为湿边温度端元。该方法的数学表达式如下：

其中LST_w代表土地覆盖类型“永久湿地”或“水体”的LST(K)；LST_{no_ice}代表没有“冰雪”的土地覆盖类型以及上述湿地和水体上的LST(k)；T_a是实时近地面气温(K)。当LST_w存在于一个局部研究区域时，LST_w的平均值被认为是T_s^min和T_c^min。如果不存在LST_w且存在LST_{no_ice}<T_a，则将LST_{no_ice}的最小值设置为T_s^min和T_c^min。对于其他情况，T_a被设为T_s^min。

本申请引入Biome-BGC模型中的平行阻抗的方式，计算干边上的温度端元。具体地，植被冠层上的空气动力阻抗r_a^c使用以下方程式确定：

式中，rh是植被表面的对流传热阻抗；rr是辐射传热阻抗；r_a^c是干燥植被冠层的空气动力学阻抗；gl_sh为叶尺度边界层导度(m/s)，其值根据土地覆盖类型确定。对于常绿针叶林、常绿阔叶林、落叶针叶林、落叶阔叶林、混交林来说，gl_sh＝0.01m/s。对于封闭灌木、开放灌木、草地、农田、贫瘠或稀疏的植被，gl_sh＝0.04m/s(Running>是干燥裸土组分的空气动力学阻抗，则使用如下方程式计算：

其中，r_hs是裸土表面的对流传热阻抗；rbl_min表示蒸汽输送的总空气动力阻力(s/m)，在本算法中其值设定为95.0s/m。

基于上述空气动力学阻抗，干边上的温度端元T_s^max和T_c^max可以通过以下公式表示：

其中α_c和α_s是植被组分和裸土组分短波的反照率(0.3-5.0μm)；ε_c和ε_s是植被组分和裸土组分的宽波段发射率(8-14μm)，在没有测量的情况下，可以从(Rubio等，2003)的查找表中获得ε_c和ε_s。为了方便起见，ε_c和ε_s的值被设为查找表中的均值，分别为0.983和0.959。r_a^c是干燥植被冠层的空气动力学阻抗，是干燥裸土组分的空气动力学阻抗。n_s是分数系数(无单位)，默认为0.35(Long和Singh>d是下行短波辐射(W/m2)。σ是Stefan-Boltzmann常数，其值约为5.67×10-8W/(m2·K4)。ρ是空气密度(kg/m3)，约为1.225kg/m3。c_p是空气比热，可近似设定为1006J/(kg·K)。ε_a是大气发射率，它可以近似地由空气温度决定，比如ε_a＝1.0-0.261exp[-7.77×10^-4(273-T_a)2](Mu，Zhao，andRunning>

上面的式4、5由Biome-BGC模型得到。

土壤湿度指数计算

土壤湿度指数SMI被定义为一种条件指数，用于测量LST/FVC空间中的任一点相对于LST/FVC空间干边和湿边的距离。具体表达式如下：

LST是遥感地表温度，单位为K，LST_dry和LST_wet是干边和湿边的温度。SMI在0到1之间变化，其中0表示极度干燥的条件，1表示极度潮湿的条件。为了评价本申请提出的方法，分别利用本申请的方法和一个理论计算方法，计算SMI，并将得到SMI与航空遥感实测的土壤湿度做对比。

下面，以下面的研究区域和数据进行计算和比较

研究区域与数据

研究区域

这项研究中收集了2012年土壤湿度主被动验证实验(SMAPVEX12)第一版本(2017)主被动L波段传感器(PALS)的土壤湿度数据。SMAPVEX12是2012年在加拿大进行的基于飞机的现场试验，它旨在满足土壤湿度算法验证和土壤湿度主动被动(SMAP)卫星任务产品的要求。从6月7日到7月19日，这次试验持续了6周。在此期间，美国国家航空航天局(NASA)每周飞行两次，两架飞机携带仪器，像SMAP卫星一样可以提供高分辨率L波段辐射计和雷达产品。在空中采集的同时，地面工作人员收集了土壤水分数据和一些其他相关参数。通过许多科学家的努力，现在可以从国家冰雪数据中心(NSIDC)获得SMAPVEX12的各种数据。有关SMAPVEX12的具体信息可在一些文献中找到(Colliander等2016年；McNairn等人2015)。利用高空飞行的PALS亮度-温度观测，获取土壤湿度数据。该算法使用传统τ-ω-h算法的变体，回归土壤湿度的准确度报告为RMSD(均方根差)为0.058>3/m³，UBRMSD(无偏RMSD)为0.056>3/m³，Pearson相关系数为0.87。

研究数据

(1)MODIS数据

本研究使用了四种modis产品，它们是：(1)MODIS/Terra和Aqua Albedo每日L3全球500米SIN Grid V006(MCD43A3)；(2)MODIS/Aqua和Terra陆地表面温度/发射率每日L3全球1公里SIN Grid V006(MYD11A1)；(3)MODIS/Terra和Aqua叶面积指数/FPAR 8天L4全球500米SIN Grid V006(MCD15A2H)；(4)MODIS/Aqua Aerosol 5-Min L2 Swath 3km(MYD04_3K)MYD11A1用于获得LST，其标准为“像素产生”，“良好数据质量”，“平均发射率误差≤0.02”和“平均LST误差≤2K”。使用MCD15A2H获得LAI(叶面积指数)并使用f_v＝1-exp(-k_parLAI)的等式计算FVC，其中f_v代表FVC，k_par是假设的消光系数，取0.5，适用于具有球形叶角分布的冠层。MCD43A3用于获得陆面夜晚和白天的反照率。在MYD04_3K的帮助下，利用这些反照率计算陆地蓝天反照率及其对土壤和植被成分的分解。蓝天反照率可以估计为白天和黑天反照率α＝f_diffα_white+(1-f_diff)α_black的线性组合，其中α是蓝天反照率；f_diff是散射分数；α_white和α_black是白天和夜晚的短波反照率，f_diff可以通过解析大气辐射传输模型结果的查找表来确定。在本次研究中，利用了6S模型建立的LUT，LUT通过太阳天顶角(0到89度，1度为单位)、光学深度(0到1.0.02度单位)和气溶胶模型类型(大陆和海洋)进行索引。从MYD04_3K产品获得太阳天顶角和气溶胶光学深度。

到达传感器的短波辐射可以假定为来自植被和土壤组分的辐射加权和。这一假设也在以往将地表反照率分离为植被和土壤组分的研究中得到了验证(Long和Singh 2012b；Zhang等。2005年；Houldcroft等2009)。基于这一假设，植被/土壤混合表面上的α可视为植被冠层反照率(α_c)和裸地反照率(α_s)的线性组合(Houldcroft等2009)即α＝f_vα_c+(1-f_v)α_s。在本次研究中，FVC高于0.9的像素上的数据平均为α_c。然后，通过求解上述反照率的线性组合来计算α_s

(2)NLDAS-2强制数据集

选择NLDAS-2强制数据集是因为它具有更高的时间(每小时)和空间分辨率(～12.5km)。该数据集来源于美国国家环境预测中心(NCEP)北美区域再分析(NARR)。NARR分析场的空间分辨率为32km，时间频率为3小时。它们在空间上插入到NLDAS 1/8度网格中，然后在时间上按照NLDAS-2小时频率进行分解。此外，地面气压和气温等气象变量也作了垂直调整，以解释NARR和NLDAS地形高度场之间的垂直差异。有关处理方法的更多详细信息，请参见(Cosgrove等2003)。.NLDAS-2的每小时陆面强迫场被分为两个GRIB文件“文件A”和“文件B”，这里只使用“文件A”，文件中需要三个变量，包括表面压力(Pa)、2米处的空气温度(K)和向下短波辐射(w/m2)。

结果分析

与其他算法比较

图2到图6是6月27日、7月5日、7月13日、7月17日四天的土壤湿度及三种不同方法计算的土壤湿度指数。在土壤湿度较低的土壤条件下，三种方法计算的SMI与SM的分布基本一致。在土壤湿度较高的一些位置，尤其是7月13日Long和Sun方法的SMI明显比本申请的方法的偏高且与SM的分布情况很不一致。Long和Sun方法的湿边高于本申请的方法的湿边，Long和Sun方法的T_s^min和T_c^min高于本申请的方法。本申请的方法的T_s^min和T_c^min相对更低，且计算得到的SMI与SM的分布具有更好的一致性。

在土壤湿度较高的一些位置Long和Sun方法的SMI明显比本申请的方法偏高是因为没有考虑土壤湿度较高时会出现湿边温度小于气温的情况。本申请的方法分条件讨论了不同土地覆盖类型T_s^min和T_c^min的不同计算方法。将研究区域分为有存在“永久湿地”或“水体”、不存在“永久湿地”或“水体”且地表温度小于气温和其他类型三种条件来分别去确定研究区域可能的湿边最低值。图6(d)直方图可以看出LST-T_a值存在大部分小于零的部分，导致图1中7月13日Long和Sun方法的SMI与SM的分布出现较大偏差。本申请的方法不仅有与Long和Sun方法可比的计算结果，并且细分了湿边温度的计算方法，使湿边温度更合理，也增加了LST/FVC空间模型的使用性。

与土壤湿度比较

图7展示了7月5日，7月13日，7月17日天的土壤湿度和各种SMI比较的散点图，横轴是土壤湿度，纵轴是基于各种算法计算出来的SMI，表1给出了三种不同方法的样本值在p＝0.05时与土壤湿度的相关系数。本申请的方法的相关系数普遍高于其他两种方法，即其分布与SM的分布的一致性强于Long和Sun方法，这也与图1中的结果相吻合。同时本申请的方法的样本值明显多于其他两种方法，体现了本申请的方法计算过程需要的参数相对较少数据易于获取的优点。

表1各个方法的SMI与SM的相关系数及样本值

其中，SMI_semi代表本申请的方法所计算的SMI，SMI_long代表long的方法所计算的SMI，SMI_sun代表sun的方法所计算的SMI。

图8展示了所有时间点上，三种不同方法的相关系数与土壤湿度平均值的变化情况。坐标横轴是时间点，左侧纵轴是平均土壤湿度与柱状图对应，右侧纵轴是相关系数与折线点对应。在土壤湿度一样的条件下，本申请的方法与Long和Sun方法相比，具有相当甚至更好的相关性。但在土壤湿度比较大(大于0.2)的条件下，三种方法的相关系数相当且都比较低。并且在土壤湿度在0.15—0.2之间相关系数接近或大于0.4，相关性比较高，且本申请的方法的最高。本申请的方法在所有日期都有研究数据，而Sun方法因为需要的参数太多导致缺少研究数据也最多。

图8显示，在干边，本申请的方法的T_c^max和T_c^max与Sun和Long方法一致，并且与土壤湿度分布高度吻合。同时在计算过程中引入Biome-BGC模型和对流(rh)、辐射(rr)传热的平行阻力、gl_sh等易于获取的参数，简化计算公式，提升了运算效率。在湿边，充分考虑了地表覆盖类型的复杂情况，增加了湿边在地表温度小于气温情况下湿边的合理性。

敏感性分析

为了更好地理解计算SMI的本申请的方法，进行了灵敏度分析。将新计算的温度终值与初始值(T_s^max和T_c^max)的相对变化与新输入变量相对于其初始值的相对变化进行了比较。输入变量包括T_a，P_a，S_d，α_s和α_c,T_a在(T_a-2K，T_a+2K)中变化，以0.5K为单位。其他参数在初始值的±20％范围内变化，每次变化2％。2012年7月5日，在一个像素(纬度＝49°46'38.17“N，经度＝98°12'36.30”W)上选择了初始场景，此处LST＝303.67K，FVC＝0.82，T_a＝28.596℃，P_a＝97896.70Pa，α_s＝0.1，ndα_c＝0.17。

图9(a)是计算出的T_s^max相对于P_a、α_s和S_d的灵敏度分析，(b)是T_c^max相对于P_a、α_c和S_d的灵敏度分析，(c)中是T_s^max和T_c^max相对于T_a的灵敏度分析。

图9显示了灵敏度分析结果，表明温度极值(即T_c^max和T_s^max)对T_a最敏感，然后是S_d。它们都与温度极值成正相关，T_a增加2K可导致温度极值增加近2K，S_d±20％的相对变化可以使这些端部成员极值发生±4K的变化。P_a对T_c^max没有影响，而与T_s^max呈正相关，P_a的相对变化为±20％时，T_s^max的相对变化约为±2.5K。两个温度极值对土壤和植被反照率不敏感。土壤反照率和植被反照率分别与T_s^max和T_c^max呈负相关,土壤反照率的±20％变化导致T_s^max值约为±0.5K，植被反照率的相同变化导致T_c^max值约为±1K扰动。

本申请的方法有几个明显的优点：

首先，它对于正常的网格化气象数据集更实用，因为它只需要三个气象变量，即气温，气压和下行短波辐射。有许多具有高空间分辨率和时间分辨率，可以供公众使用的网格化气象数据集(Walton和Hall 2018)，比如例如NLDAS-2强迫数据集，Prism(独立斜坡模型上的参数高程回归)(Daly等2008年)、Gridmet(网格地面气象数据)(Abatzoglou 2013年)和WRF(天气研究和预测)(Walton等2017)等。将这些网格化气象数据集与遥感LST/FVC相结合，可以改善半经验方法的土壤水分监测或微波SM降尺度。

其次，本申请的方法是一种非迭代方法，这可以避免迭代算法的一些缺点，节省计算时间，避免不收敛，提高计算效率。

第三，本申请的方法开发了一种“逐步的方法”来确定湿边。这样可以避免气温高于地表温度的现象，SMI_long和SMI_sun在这种情况下失效。

最后，利用Biome-BGC模型避免了非常复杂的空气动力学阻抗计算。