技术领域
[0001] 本
发明实施例涉及生态
水文学与遥感领域，具体涉及一种流域植被干旱发生遥感预警与缺水估算方法。
背景技术
[0002] 干旱区指降水稀少、
蒸发量偏大、产流量很少的地区。植被
生态系统水分胁迫通常指
土壤缺水而明显抑制植被生长的现象。淹水、
冰冻、高温或盐渍等也能引起水分胁迫。干旱缺水引起的水分胁迫是最常见的，也是对
植物产水量影响最大的。目前，水资源协调和管理多是以流域、子流域、二级子流域为单元开展的。
[0003] 面向不同的流域尺度的植被生态水分胁迫遥感识别的模型与方法，现有的植被生态系统水分胁迫遥感发现的方法主要包括
微波或热红外反演土壤水分、蒸散发、降水分析法、植被生理分析法等。其中，微波或热红外反演土壤水分探测土壤水分的深度有限、受遥感影像的空间
分辨率影响、受反演结果不稳定影响等；蒸散发遥感反演的方法受限的原因是得到的蒸散发量是卫星过境的瞬时值，并在此
基础上进行推算累加得到的，因此受卫星观测
频率、反演的模型参量、遥感影像的空间分辨率、天气变化(日照、降水、
云量等)等因素影响，具有较大的不确定性；降水分析方法的主要问题是主要从大区域和气象的
角度，在一定时间段内的降水量和分布进行分析，对于不同尺度的流域或子流域缺乏考虑降水产生的径流、蒸散发等重要因素；植被生理分析法主要通过遥感方法识别植被缺水导致的叶面积指数、冠层
含水量、叶表面
温度等生理参量变化，该方法的问题反演结果的
稳定性和
精度都有待提高。
[0004] 上述方法多适用于相对均一的下垫面情况。同时，现有的植被生态水分胁迫遥感识别的方法，多是以
像素为单元的，而水资源协调与管理多是以不同尺度的流域为单元开展的。
[0005] 故存在，流域不同尺度的气象干旱预警中，气象干旱预警方法由于缺乏考虑土壤水分变化，及灌区和非灌区作物、林、草等植被类型及分布等因素导致精度不高，而且不能有效的识别流域尺度干旱的差异；干旱导致的缺水量估算方法中，采用单一的降水分析方法不能有效识别灌区、非灌区、林草等耗水差异存在明显的
缺陷，但目前的遥感蒸散发的方法主要通过卫星反演的瞬时值进行时间尺度扩展，存在较大的不确定性，从而导致多尺度流域缺水估算精度不高。
发明内容
[0006] 为此，本发明实施例提供一种流域植被干旱发生遥感预警与缺水估算方法，以解决现有流域植被干旱不能进行胁迫预警，多尺度流域缺水估算精度不高的问题。
[0007] 为了实现上述目的，本发明实施例提供如下技术方案：
[0008] 本发明实施例公开了一种流域植被干旱发生遥感预警与缺水估算方法，所述方法包括：
[0009] 步骤1，空间
数据处理与信息提取，提取多尺度流域信息、植被类型信息、径流信息、降水信息和土壤水分信息；
[0010] 步骤2，进行干旱胁迫预警，计算水资源量变化量和同期植被变化量，根据瞬时蒸散发指数进行干旱预警；
[0011] 步骤3，对植被干旱缺水量进行估算，根据特定时间段内历史多年平均耗水量估算、监测年特定时间段耗水量估算，进行监测特定时间段缺水量估算。
[0012] 进一步地，所述步骤1中进行空间数据处理与信息提取，包括：多尺度流域提取与地面数据获取模
块、潜在蒸发量计算模块、实际蒸发量计算模块、植被信息提取模块和生态系统耗水量计算模块；
[0013] 所述多尺度流域提取与地面数据获取模块，通过选择数字高程模型DEM数据，提取n级流域地面数据、n-1级流域地面数据……1级流域地面数据，其中1级流域指全流域，2、3….n级流域为不同级别的子流域，地面数据包括径流量和土壤湿度；
[0014] 所述潜在蒸发量计算模块，利用第一预设方法计算植物系数Kc，Kc的计算公式为：
[0015] Kc＝Kcb+0.05
[0016] Kcb＝Kcbmin+(Kcbfull-Kcbmin)×min(1,2×fcell1/(1+h))
[0017] fcell＝fc/sinη
[0018]
[0019]
[0020] Kcb,h＝1.0+0.1×h h≤2m
[0021] Kcb,h＝1.2 h＞2m
[0022] 式中，Kcbmin为裸地最小植物系数，Kcbfull为全
覆盖条件下且叶面指数低于3时的基本植物系数，fcell为有效植被覆盖度，fc遥感
图像空间离散化反演所得为植被覆盖率，η为太阳高于海平面的角度，δ1为太阳赤纬， 为从DEM图像提取的太阳纬度，Kcb,h为标准湿润和
风速下的全覆盖植被的基本植物系数，U2为植物生长期内2m高
位置处的平均风速，RHmin为植物生长器的平均最小
相对湿度，h为平均植被最大高度；
[0023] 利用第二预设方法通过以下公式求取植物系数Kc所表示的区域内的潜在
蒸散量ETp0，
[0024]
[0025]
[0026]
[0027]
[0028] ed＝RH×es
[0029]
[0030]
[0031] Rn＝(1-αshort)Rswd↓+Rlwd↓-Rlwd↑-(1-ε0)Rlwd↓
[0032] ＝(1-αshort)Q+ε0Rlwd↓-ε0σTs4
[0033] Q＝Q0(a+bS1)ed＝RH×es
[0034] 式中，Δ为气温Ta时的饱和水汽压曲线斜率，Ta为月平均温度，γ为干湿常数，Cp为空气定压
比热，Pr为
大气压，ε为水汽分子量与干空气分子量之比，λ为蒸发
潜热，H为根据DEM图获取的海拔高度，eS为平均饱和水汽压，ed为实际水汽压，RH为空气相对湿度，Uz为z高度处观测到的风速,G为地表热通量，Ts为地表温度，α为地表反照率，NDVI为通过遥感图像反演所得的归一化植被指数，Q为太阳总
辐射量，Q0为太阳总辐射量初始值，a与b均为经验系数，S1为同期的日照百分率，Rn为净辐射量，αshort为短波反照率，Rswd↓为大气下行短波辐射，Rlwd↓为大气向下的长波辐射，Rlwd↑为大气向上的长波辐射，ε0为地表宽波段比辐射率，σ的为Stefan-Boltzmann常数，值为4.903×10-9MJ m-2K-4。
[0035] 进一步地，所述实际蒸发量计算模块，利用第三预设方法求取土壤蒸发量Ea，按照以下公式估算出绿水流估算值Z，
[0036] ETp＝Kc·ETp0
[0037]
[0038] Z＝Ecan,a+Ta+Ea
[0039] 式中，ETp为植物潜在蒸散量，Ecan,a为植物冠层截留蒸量，cp为通过遥感图像空间离散化像元反演得到植物覆盖率，δ为湿润叶面面积率，Ta为植物蒸腾量，Si为i层根系吸水量，n为植物根系总层数。
[0040] 进一步地，所述植被信息提取模块对灌区和非灌区进行区分，确定复种与种植结构信息。
[0041] 进一步地，所述确定复种与种植结构信息包括以下步骤：
[0042] 判定作物种植与成熟的时间点，通过历史数据，选择降水和土壤水分正常的年份，建立时间序列的归一化植被指数NDVI曲线，判定植物种植与成熟的时间
节点；
[0043] 建立同像素的多期作物生长曲线，提取作物种植结构和复种信息，构建不同类型的作物从种植到收割的生育周期NDVI变化曲线，t-NDVI曲线，重构时间序列的NDVI数据，建立时间t-NDVI曲线，根据t-NDVI曲线，在时间序列的NDVI数据中选择样本，采取监督分类的方法，获得作物种植结构信息、不同作物的复种信息。
[0044] 进一步地，所述灌区和非灌区的区分方法为：判断出湿润期和干旱期，通过特定时间段内土壤水分变化分别和降水正常年、降水湿润年以及降水干旱期的NDVI均值大小，判断灌区或非灌区。
[0045] 进一步地，所述生态系统耗水量计算模块通过以下公式计算t2-t1时间段内植被生态耗水量，
[0046] ECt2-t1＝ΔRt2-t1-ΔWt2-t1+ΔSt2-t1+ΔPt2-t1
[0047]
[0048] SWi＝(RSWIR-RG)/(RSWIR-RG)≥Δb1
[0049] 式中ECt2-t1为t2-t1时间段内植被生态耗水量，ΔRt2-t1为t1-t2时间段内的流域产生的径流量，采用地面测量数据；ΔWt2-t1为t1-t2时间段内流域水面蒸发量，ΔSt2-t1为t1-t2时间段内流域土壤水分变化量，ΔPt2-t1为t1-t2时间段内流域降水量；EWi为时段i的水面蒸发量，采用地面测量数据；SWi为时段i的水面面积，由遥感数据获得；RSWIR为遥感的短波红外波段反射率，RG为绿波段的反射率；Δb1为预设
阈值。
[0050] 进一步地，所述步骤2中，干旱胁迫预警的方法为：计算水资源变化指标、同期植被指数变化指标和历史同期指标，计算瞬时蒸散发指数，在统计的时间段内将一定深度的土壤水资源变化量与预设阈值进行比较，判断是否进行干旱预警。
[0051] 进一步地，所述步骤3中，植被干旱缺水量估算方法为：估算历史多年特定时间段生态平均耗水量，估算监测年特定时间段的生态耗水量，监测特定时间段的缺水量为历史多年特定时间段生态平均耗水量与监测年特定时间段的生态耗水量的差值。
[0052] 本发明实施例具有如下优点：
[0053] 本发明实施例公开了一种流域植被干旱发生遥感预警与缺水估算方法，对空间数据进行处理并进行信息提取，获取多尺度流域信息、植被类型信息、径流信息、降水信息和土壤水分信息，对地形进行分析，充分考虑流域的尺度性和整体性以及不同植被空间耗水差异，在统计的时间段内将一定深度的土壤水资源变化量与预设阈值进行比较，进行干旱预警，根据历史多年特定时间段生态平均耗水量与监测年特定时间段的生态耗水量估算植被干旱缺水量，考虑气象数据和地面观测数据的有效性，保证较高精度的多尺度流域生态系统干旱预警与缺水量估算，便于及时对干旱区域进行浇灌，缓解干旱状态。
附图说明
[0054] 为了更清楚地说明本发明的实施方式或
现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。
[0055] 本
说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。
[0056] 图1为本发明实施例提供的一种流域植被干旱发生遥感预警与缺水估算方法
流程图；
具体实施方式
[0057] 以下由特定的具体实施例说明本发明的实施方式，熟悉此技术的人士可由本说明书所揭露的内容轻易地了解本发明的其他优点及功效，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0058] 实施例
[0059] 本实施例公开了一种流域植被干旱发生遥感预警与缺水估算方法，包括：
[0060] 步骤1，空间数据处理与信息提取，提取多尺度流域信息、植被类型信息、径流信息、降水信息和土壤水分信息；
[0061] 步骤2，进行干旱胁迫预警，计算水资源量变化量和同期植被变化量，根据瞬时蒸散发指数进行干旱预警；
[0062] 步骤3，对植被干旱缺水量进行估算，根据特定时间段内历史多年平均耗水量估算、监测年特定时间段耗水量估算，进行监测特定时间段缺水量估算。
[0063] 步骤1中进行空间数据处理与信息提取，包括：多尺度流域提取与地面数据获取模块、潜在蒸发量计算模块、实际蒸发量计算模块、植被信息提取模块和生态系统耗水量计算模块；
[0064] 所述多尺度流域提取与地面数据获取模块，通过选择数字高程模型DEM数据，提取n级流域地面数据、n-1级流域地面数据……1级流域地面数据，其中1级流域指全流域，2、3….n级流域为不同级别的子流域，地面数据包括径流量和土壤湿度；
[0065] 所述潜在蒸发量计算模块，利用第一预设方法计算植物系数Kc，Kc的计算公式为：
[0066] Kc＝Kcb+0.05
[0067] Kcb＝Kcbmin+(Kcbfull-Kcbmin)×min(1,2×fcell1/(1+h))
[0068] fcell＝fc/sinη
[0069]
[0070]
[0071] Kcb,h＝1.0+0.1×h h≤2m
[0072] Kcb,h＝1.2 h＞2m
[0073] 式中，Kcbmin为裸地最小植物系数，Kcbfull为全覆盖条件下且叶面指数低于3时的基本植物系数，fcell为有效植被覆盖度，fc遥感图像空间离散化反演所得为植被覆盖率，η为太阳高于海平面的角度，δ1为太阳赤纬， 为从DEM图像提取的太阳纬度，Kcb,h为标准湿润和风速下的全覆盖植被的基本植物系数，U2为植物生长期内2m高位置处的平均风速，RHmin为植物生长器的平均最小相对湿度，h为平均植被最大高度；
[0074] 利用第二预设方法通过以下公式求取植物系数Kc所表示的区域内的潜在蒸散量ETp0，
[0075]
[0076]
[0077]
[0078]
[0079] ed＝RH×es
[0080]
[0081]
[0082] Rn＝(1-αshort)Rswd↓+Rlwd↓-Rlwd↑-(1-ε0)Rlwd↓
[0083] ＝(1-αshort)Q+ε0Rlwd↓-ε0σTs4
[0084] Q＝Q0(a+bS1)ed＝RH×es
[0085] 式中，Δ为气温Ta时的饱和水汽压曲线斜率，Ta为月平均温度，γ为干湿常数，Cp为空气定压比热，Pr为大气压，ε为水汽分子量与干空气分子量之比，λ为蒸发潜热，H为根据DEM图获取的海拔高度，eS为平均饱和水汽压，ed为实际水汽压，RH为空气相对湿度，Uz为z高度处观测到的风速,G为地表热通量，Ts为地表温度，α为地表反照率，NDVI为通过遥感图像反演所得的归一化植被指数，Q为太阳总辐射量，Q0为太阳总辐射量初始值，a与b均为经验系数，S1为同期的日照百分率，Rn为净辐射量，αshort为短波反照率，Rswd↓为大气下行短波辐射，Rlwd↓为大气向下的长波辐射，Rlwd↑为大气向上的长波辐射，ε0为地表宽波段比辐射率，σ的为Stefan-Boltzmann常数，值为4.903×10-9MJ m-2K-4。
[0086] 实际蒸发量计算模块，利用第三预设方法求取土壤蒸发量Ea，按照以下公式估算出绿水流估算值Z，
[0087] ETp＝Kc·ETp0
[0088]
[0089] Z＝Ecan,a+Ta+Ea
[0090] 式中，ETp为植物潜在蒸散量，Ecan,a为植物冠层截留蒸量，cp为通过遥感图像空间离散化像元反演得到植物覆盖率，δ为湿润叶面面积率，Ta为植物蒸腾量，Si为i层根系吸水量，n为植物根系总层数。
[0091] 植被信息提取模块对灌区和非灌区进行区分，确定复种与种植结构信息，所述确定复种与种植结构信息包括以下步骤：
[0092] 判定作物种植与成熟的时间点，通过历史数据，选择降水和土壤水分正常的年份，建立时间序列的归一化植被指数NDVI曲线，判定植物种植与成熟的时间节点；
[0093] 在历史数据中， 则判定为水分供给正常年，选择该年的时间序列NDVI曲线数据集。式中，ΔPSj为时段j内的土壤和水分变化量；ΔPyi为时间段j内的降水量，ΔSyi为时间段j内的土壤水分减少量，如果呈增加趋势，则ΔSyi＝0，n为时段j内覆盖同一区域的NDVI获取的次数。
[0094] 判断种植和成熟的时间节点，当NDVIt-1≤NDVIt＜NDVIt+1或同时满足NDVIt＜NDVIt+1和NDVIt≤NDVIt-1，且NDVIt+1-NDVIt≥Δa1和a2≥NDVIt≥a1，时，则判定t时间点为种植的起始点，式中NDVIt、NDVIt-1和NDVIt+1分别为时间点t和t-1和t+1的NDVI值；Δa1、a1、a2分别为预设阈值；
[0095] 当NDVIt+1≤NDVIt＜NDVIt-1时，且a4≥NDVIt≥a3，NDVIt-1-NDVIt≥Δa2，NDVIt-NDVIt+1≥Δa3则判定t时间点为开始成熟期。式中，a3、a4、Δa2、Δa3分别为预设阈值。
[0096] 建立同像素的多期作物生长曲线，提取作物种植结构和复种信息，构建不同类型的作物从种植到收割的生育周期NDVI变化曲线，t-NDVI曲线，重构时间序列的NDVI数据，建立时间t-NDVI曲线，根据t-NDVI曲线，在时间序列的NDVI数据中选择样本，采取监督分类的方法，获得作物种植结构信息、不同作物的复种信息。
[0097] 灌区和非灌区的区分方法为：进行湿润期和干旱期的选择，且 时，判定为湿润期，式中，
为多年历史数据，在时段i的降水和土壤水分变化量的均值，ΔPyij为第j年时段i内的降水量，ΔSyij为第j年时段i内的降水量，n为统计的历史数据的年数， 为时段i土壤水分变化和降水正常年的NDVI均值， 为时段i土壤水分变化和降水湿润期的
NDVI均值，Δa5为预设值；
当 且 则判断为干
旱期，式中 为时段i土壤水分变化和降水干旱期的NDVI均值，Δa6为预设值。
[0098] 灌区与非灌区判断规则为， 则判断为灌区，否则为非灌区，式中a8和a9为预设阈值。
[0099] 生态系统耗水量计算模块通过以下公式计算t2-t1时间段内植被生态耗水量，[0100] ECt2-t1＝ΔRt2-t1-ΔWt2-t1+ΔSt2-t1+ΔPt2-t1
[0101]
[0102] SWi＝(RSWIR-RG)/(RSWIR-RG)≥Δb1
[0103] 式中ECt2-t1为t2-t1时间段内植被生态耗水量，ΔRt2-t1为t1-t2时间段内的流域产生的径流量，采用地面测量数据；ΔWt2-t1为t1-t2时间段内流域水面蒸发量，ΔSt2-t1为t1-t2时间段内流域土壤水分变化量，ΔPt2-t1为t1-t2时间段内流域降水量；EWi为时段i的水面蒸发量，采用地面测量数据；SWi为时段i的水面面积，由遥感数据获得；RSWIR为遥感的短波红外波段反射率，RG为绿波段的反射率；Δb1为预设阈值。
[0104] 步骤2中，干旱胁迫预警的方法为：计算水资源量变化指标，ΔSyi＝Syt2-Syt1式中， 为统计的时间段内一定深度的土壤水资源和降水资源变化量；i为统计的历史年份，n为统计的历史年数；ΔSyi为统计的时间段内一定深度的土壤水资源变化量，Syt2为时间节点t2一定深度的土壤水资源量，Syt1为时间节点t1一定深度的土壤水资源量；ΔPyi为某一时间段内的降水总量。
[0105] 计算同期植被指数变化指标，选择包括有林地、草地、农田、
水体、果园等遥感解译的土地利用图，选择过去斑块面积变化较小的区域。斑块选择的条件是，历史同期和监测期内斑块面积变化不大，将栅格数据矢量化，按照斑块进行如下计算；
[0106] 同一生长期内：ΔNDVIn＝NDVIPLyt2-NDVIPLyt1，ΔNDVIn为生长期(生长点到成熟点)内t1-t2时间段内的NDVI变化量，NDVIPLyt1为时间点t1的NDVI值，NDVIPLyt2为时间点t2的NDVI值。
[0107] 计算历史同期指标， 为历史同期时间区间t1-t2的NDVI变化值； 为历史同期t2时间点的NDVI平均值； 为历史同期
t1时间点的NDVI平均值。
[0108] 计算瞬时蒸散发指数，EI＝(ETsp-ETsa)/ETsp，EI为瞬时蒸散发指数；ETsp为瞬时潜在蒸散发；ETsa为瞬时实际蒸散发。
[0109] 进行干旱预警， 且ΔPSi≤Δb2、ΔNDVIn≤Δb3、EI≤Δb5，WRi＝ΔPSi/b1+ΔNDVIn/b2+EI/b3，当WRi≤c0时，进行干旱预警。式中，b2、b3、b5、c0为预设阈值。
[0110] 当DRwjk≥c2时，为流域干旱预警，DRwjk为k级流域干旱胁迫指数，Sdi为受干旱胁迫的斑块i的面积，Sk为k级流域的面积，Ska为灌区的面积。
[0111] 步骤3中，植被干旱缺水量估算方法为：估算历史多年i时段平均耗水量，式中 为历史多年i时段生态平均耗水量；ΔRoij为历史多年j年i时段径流产生量；ΔWoij为j年i时段自由水面蒸发量；ΔSoij为j年i时段土壤水分变化量；ΔPoij为j年i时段的降水总量；ΔIRoij为j年i时段
灌溉取
地下水量；j为统计的年，n为年数；
[0112] 估算监测年i时段耗水量，ECi＝ΔRi-ΔWi+ΔSi+ΔPi+ΔIRi，式中ECi为监测年i时段的生态耗水量，ΔRi为i时段的径流量，ΔWi为i时段的自由水面蒸发量，ΔSi为i时段土壤水资源变化量，ΔPi为i时段的降水量，ΔIRi为灌溉取地下水量；
[0113] 估算监测时段i的缺水量， 式中ECs为监测时段i的流域缺水量。
[0114] 本发明实施例公开的一种流域植被干旱发生遥感预警与缺水估算方法，对空间数据进行处理并进行信息提取，获取多尺度流域信息、植被类型信息、径流信息、降水信息和土壤水分信息，对地形进行分析，充分考虑流域的尺度性和整体性以及不同植被空间耗水差异，在统计的时间段内将一定深度的土壤水资源变化量与预设阈值进行比较，进行干旱预警，根据历史多年特定时间段生态平均耗水量与监测年特定时间段的生态耗水量估算植被干旱缺水量，考虑气象数据和地面观测数据的有效性，保证较高精度的多尺度流域生态系统干旱预警与缺水量估算，便于及时对干旱区域进行浇灌，缓解干旱状态。
[0115] 虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些
修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。
