CN103645295B - 一种多层土壤水分模拟方法和系统 - Google Patents

Abstract

本发明提供一种多层土壤水分模拟方法和系统，包括：构建土壤水分层均衡模型，结合遥感技术调整模型结构；利用遥感技术获取模型参数，构建流域水文空间信息数据库；利用调整后的土壤水分层均衡模型和所述流域水文空间信息数据库，开展多层土壤水分过程数值模拟。本发明能够将遥感反演方法和传统的典型土壤水分层均衡模型相耦合，构建出物理机制明确、结构简单的模型，并能保持一定的土壤水模拟精度，实现多层土壤水分模拟的空间化和动态化。

Description

一种多层土壤水分模拟方法和系统

技术领域

本发明涉及生态水文科学领域，尤其涉及一种多层土壤水分模拟方法和系统。

背景技术

土壤水分为水循环的重要组成部分，是水文学、气象学以及农业科学研究领域中的一个重要指标参数。土壤水的含量和运动变化，深刻地影响土壤内部物质的转化过程，并通过蒸发蒸腾的方式在农业生产、生态环境建设和保护方面发挥着重要作用。在全球土壤湿度计划(Global Soil Wetness Porject，GSWP)中，土壤水已被作为专门的议题来进行探讨和研究，且已经成为全球变化研究的重要焦点之一。因此，深入研究土壤水分变化机理，弄清土壤水分含量和运移规律，对加强区域水资源管理，合理利用降水资源，缓解水资源匮乏具有重要的现实意义。

所谓土壤水是指在一个大气压下，在105℃条件下能从土壤中分离出来的水分。土壤水是植物生长和生存的物质基础，它不仅影响林木、大田作物、蔬菜、果树的产量，还影响陆地表面植物的分布。

下垫面是指地球表面，包括海洋、陆地、陆地上的高原、山地、平原、森林、草原以及城市等。下垫面各部分温度、水分以及表面形状等参数均有较大差异，从而导致下垫面具有非均一性。

传统的土壤水分监测方法有经验统计法(从气象和地形因素通过经验公式和数理统计方法计算土壤水分)、应用微气象学方法(能量平衡、空气动力能量平衡和空气动力涡度相关等)、土壤水运动规律分析法(土壤水量平衡计算、零通量面、测渗学法和土壤水动力学等)。而土壤水分的遥感反演方法主要有热惯量法、温度植被指数法、微波遥感法等。

但是，现有技术中的土壤水分监测方法虽然可以模拟厚层土壤水分变化规律，但是很难满足区域尺度土壤水分连续监测；遥感反演方法虽然具有快速、客观的优势，在一定程度上解决了表层和浅层土壤湿度的反演问题，但是大都在特定的时间和地域建立的，很难满足深层土壤湿度反演的需要。因此，传统的土壤水分监测方法和遥感反演方法都具有一定的局限性，无法满足多层土壤水分模拟的空间化和动态化需求。

发明内容

(一)要解决的技术问题

本发明提供一种多层土壤水分模拟方法和系统，以解决现有技术中无法满足多层土壤水分模拟的空间化和动态化需求的技术问题。

(二)技术方案

为解决上述技术问题，本发明提供一种多层土壤水分模拟方法，其特征在于，包括：

构建土壤水分层均衡模型，结合遥感技术调整模型结构；

利用遥感技术获取模型参数，构建流域水文空间信息数据库；

利用调整后的土壤水分层均衡模型和所述流域水文空间信息数据库，开展多层土壤水分过程数值模拟。

进一步地，所述土壤水分层均衡模型包括：土壤水分平衡模型、植被截流模型、地表径流模型、蒸散发模型、根系吸水模型、层间土壤水分增量模型、根系层补给量模型中的一个或多个。

进一步地，所述构建土壤水分层均衡模型，结合遥感技术调整模型结构包括：

构建并调整土壤水分平衡模型：

降水条件下，表层、有棵间蒸发的土壤水分平衡方程为：

$W_1^{j+1}=W_1^j+{\mathrm{\ΔD}}_1^j-S_1^j,$

其余为：

非降水条件下，表层土壤水分平衡方程为：

$W_1^{j+1}=W_1^j+s\_uppersop\_{\mathrm{loss}}_2-{\mathrm{Ea}}^j-S_1^j,$

其余为：

Δs_upperrop_loss_i＝s_upperrop_loss_i+1-s_upperrop_loss_i-1(i≥2)；

其中i为土壤层次，j为时间节点，为第i层时段初根系层土壤贮水量，为第i层时段末根系层土壤贮水量，E_a为表层棵间蒸发量，ΔD_i为降水或灌水后第i层土壤水分增量，S_i为第i层作物根系的吸水量，s_upperrop_loss_i为第i层净水分上移量；

和/或，构建并调整植被截流模型：

构建林冠截留量计算方程为：

$S_v=c_v\×S_{max}\×\[1-e^{-\mathrm{\η}\frac{P}{S_{max}}}\]$

式中，S_v为累计截留量，单位mm；c_v为植被盖度，单位％，反映覆被空间分布情况；P为累积降水量，单位mm；S_max为树冠蓄水能力，单位mm，即林冠最大截留量；η为校正系数；

其中林冠最大截留量为：

S_max＝0.935+0.498×LAI-0.00575×LAI²，

校正系数：η＝0.046×LAI，LAI为叶面积指数；

和/或，构建并调整地表径流模型：

$Rs=g_1{\left(\frac{{\mathrm{AW}}_u}{{\mathrm{WM}}_u\·C}\right)}^{g_2}\·P^{\′}$

式中，Rs为地表产流量，单位mm；AW_u为表层土壤湿度，单位％；WM_u为表层土壤饱和含水量，单位％；P'为有效降雨量，单位mm，是降雨量扣除林冠截留后到达地面的净雨量；g₁与g₂是时变增益因子的有关参数，0<g₁<1，1<g₂，其中g₁为土壤饱和后径流系数，g₂为土壤水影响系数；C为覆被影响参数；

和/或，构建并调整蒸散发模型：

植被覆盖区实际蒸散量由土表实际棵间蒸发量E_a和实际蒸腾量T_a组成：

E_a＝K_s×E_p，

T_a＝K_s×T_p，

其中土壤水分胁迫系数K_s为：

K_s＝ln(A_v+1)/ln(101)，

A_v＝[(W-W_a)/(W_f-W_a)]×100，

式中，W_f为田间持水量，单位mm；W_a为风干土含水量，单位mm；A_v为相对有效含水量。

土表潜在棵间蒸发量E_p为：

LAI为叶面积指数，

植被潜在蒸腾量T_p为：

T_p＝ET_p-E_p；

潜在蒸散ET_p的表达式为：

ET_p＝K_c·ET_p0，

式中，K_c为作物系数；ET_p0为参考作物蒸散量，单位mm；

和/或，构建并调整根系吸水模型：

根系吸水量的表达式为：

$S=T_p\&times;\frac{e^{-AROOT\&times;z_2}-e^{-AROOT\&times;z_1}}{e^{-AROOT\&times;LR}-1},$

其中T_p为植被潜在蒸腾量，AROOT为描述根主要分布的参数，z₁、z₂分别为所求土壤层垂直方向上的两端坐标，单位m，LR为根系深度，单位m；

和/或，构建并调整层间土壤水分增量模型：

当进入第i层的剩余入渗量大于AD_max,i时，降水或灌溉导致的层间土壤水分增量：AD_i＝AD_max,i，其中AD_max,i为各层的最大增量，AD_max,i＝(W_f,i-W_i)·L_i，i为土壤层数，L_i为土层厚度，取10mm；

否则，层间土壤水分增量等于剩余入渗水量：

$Gd=\left\{\begin{array}{c}F-\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i(F>\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i)\\ 0\end{array}\right.,$

其中，F＝P+I-S_V-Rs，P为期内降水量，单位mm，I为灌水量，在非灌区及天然植被区为0，单位mm；S_V为期内累计植被截留量，单位mm；Rs为期内地表径流量，单位mm；

和/或，构建并调整根系层补给量模型：

因上层根系吸水损失的水分表示为：

S_uppersop_loss_i＝K_si*ET_p*(number_layer-i)/depth_sum

其中，number_layer＝10，depth_sum＝100,单位为mm。

进一步地，

所述参考作物蒸散量ET_p0的计算方法为：

${\mathrm{ET}}_{p0}=\frac{0.408\mathrm{\&Delta;}(R_n-G)+\mathrm{\&gamma;}\frac{900}{T_a+273}{U}_2(e_S-e_d)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+0.34U_2)},$

式中：△是气温T时的饱和水汽压曲线斜率，单位kPa·℃^-1：

T是日平均温度，单位℃；

γ是干湿表常数，单位kPa·℃^-1：C_p是空气定压比热，取1.013×10^-3MJ·kg^-1·℃^-1；P_r是大气压，单位kPa；λ是蒸发潜热，取2.45MJ·kg^-1；ε是水汽分子量与干空气分子量之比，取0.622；

e_S是气温T下的饱和水汽压，单位kPa：

e_d是实际水汽压，单位kPa：e_d＝RH×e_s，RH是空气相对湿度，单位％；

U₂是2米高度处风速，单位m·s^-1：U_z为z高度处观测到的风速，单位m·s^-1；

R_n是净辐射，单位MJ·m^-2·d^-1：

$R_n=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_n\left(t\right)dt=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_{n\_max}sin\&lsqb;\left(\frac{t-t_{rise}}{t_{set}-t_{rise}}\right)\mathrm{\&pi;}\&rsqb;dt,$

t_rise，t_set分别为日出时间和日落时间，R_{n_max}为日最大净辐射，t为卫星过境时刻的当地时间；

G是土壤热通量，单位MJ·m^-2·d^-1：有植被覆盖的地面为：

G＝R_n[Γ_c+(1-f_c)(Γ_s-Γ_c)]

式中：全植被覆盖下，土壤热通量与净辐射的比值Γ_c＝0.05；裸地情况下，土壤热通量与净辐射比值Γ_s＝0.315，f_c为植被覆盖率；

对于水体和冰雪，G＝0.5R_n；

所述作物系数K_c的计算方法为：

非生长期的K_c小于0.4；

生长期的K_c为：K_c＝K_cb+0.05，K_cb为非完全覆盖条件下的基本植物系数：

式中，K_cbmin为在缺少植被情况下的裸地最小植物系数，取0.15-0.20；K_cbfull为全覆盖条件下，叶面积指数大于3的基本植物系数；h为植被高度，单位m；f_cell为有效植被覆盖度。

进一步地，

所述模型参数包括：气象参数、下垫面植被参数、下垫面土壤参数、经验参数和其它参数中的一种或几种；

所述气象参数包括：降雨量、灌溉量、气温、地表温度、2米处风速、大气压、空气相对湿度、生长期内平均风速和平均最小相对湿度中的一种或几种；

所述下垫面植被参数包括：叶面积指数、植被盖度、地表反照率、波段31发射率、波段32发射率、植被高度和平均植被最大高度中的一种或几种；

所述下垫面土壤参数包括：初始土壤含水量、土壤田间含水量、土壤饱和含水量、土壤萎蔫含水量中的一种或几种；

所述经验参数包括：地表径流参数1、地表径流参数2和地表覆被因子中的一种或几种；

所述其它参数包括：中分辨率成像光谱仪过境时间、日出时间、日落时间和地理纬度中的一种或几种。

另一方面，本发明还提供一种多层土壤水分模拟系统，包括：模型构建单元、遥感参数单元和模拟单元，模型构建单元和遥感参数单元分别与模拟单元相连，其中：

模型构建单元，用于构建土壤水分层均衡模型，结合遥感技术调整模型结构；

遥感参数单元，用于利用遥感技术获取模型参数，构建流域水文空间信息数据库；

模拟单元，用于利用调整后的土壤水分层均衡模型和所述流域水文空间信息数据库，开展多层土壤水分过程数值模拟。

进一步地，所述模型构建单元包括：

模型分类子单元，用于将土壤水分层均衡模型分为：土壤水分平衡模型、植被截流模型、地表径流模型、蒸散发模型、根系吸水模型、层间土壤水分增量模型、根系层补给量模型中的一种或多种类型。

进一步地，所述模型构建单元包括：

土壤水分平衡模型子单元，用于构建并调整土壤水分平衡模型：

降水条件下，表层、有棵间蒸发的土壤水分平衡方程为：

$W_1^{j+1}=W_1^j+{\mathrm{\&Delta;D}}_1^j-S_1^j,$

其余为：

非降水条件下，表层土壤水分平衡方程为：

$W_1^{j+1}=W_1^j+s\_uppersop\_{\mathrm{loss}}_2-{\mathrm{Ea}}^j-S_1^j,$

其余为：

Δs_upperrop_loss_i＝s_upperrop_loss_i+1-s_upperrop_loss_i-1(i≥2)；

其中i为土壤层次，j为时间节点，为第i层时段初根系层土壤贮水量，为第i层时段末根系层土壤贮水量，E_a为表层棵间蒸发量，ΔD_i为降水或灌水后第i层土壤水分增量，S_i为第i层作物根系的吸水量，s_upperrop_loss_i为第i层净水分上移量；

和/或，植被截流模型子单元，用于构建并调整植被截流模型：

构建林冠截留量计算方程为：

$S_v=c_v\&times;S_{max}\&times;\&lsqb;1-e^{-\mathrm{\&eta;}\frac{P}{S_{max}}}\&rsqb;$

式中，S_v为累计截留量，单位mm；c_v为植被盖度，单位％，反映覆被空间分布情况；P为累积降水量，单位mm；S_max为树冠蓄水能力，单位mm，即林冠最大截留量；η为校正系数；

其中林冠最大截留量为：

S_max＝0.935+0.498×LAI-0.00575×LAI²，

校正系数：η＝0.046×LAI，LAI为叶面积指数；

和/或，地表径流模型子单元，用于构建并调整地表径流模型：

$Rs=g_1{\left(\frac{{\mathrm{AW}}_u}{{\mathrm{WM}}_u\&CenterDot;C}\right)}^{g_2}\&CenterDot;P^{\&prime;}$

式中，Rs为地表产流量，单位mm；AW_u为表层土壤湿度，单位％；WM_u为表层土壤饱和含水量，单位％；P'为有效降雨量，单位mm，是降雨量扣除林冠截留后到达地面的净雨量；g₁与g₂是时变增益因子的有关参数，0<g₁<1，1<g₂，其中g₁为土壤饱和后径流系数，g₂为土壤水影响系数；C为覆被影响参数；

和/或，蒸散发模型子单元，用于构建并调整蒸散发模型：

植被覆盖区实际蒸散量由土表实际棵间蒸发量E_a和实际蒸腾量T_a组成，公式为：

E_a＝K_s×E_p，

T_a＝K_s×T_p，

其中土壤水分胁迫系数K_s为：

K_s＝ln(A_v+1)/ln(101)，

A_v＝[(W-W_a)/(W_f-W_a)]×100，

式中，W_f为田间持水量，单位mm；W_a为风干土含水量，单位mm；A_v为相对有效含水量。

土表潜在棵间蒸发量E_p为：

LAI为叶面积指数，

植被潜在蒸腾量T_p为：

T_p＝ET_p-E_p；

潜在蒸散ET_p的表达式为：

ET_p＝K_c·ET_p0，

式中，K_c为作物系数；ET_p0为参考作物蒸散量，单位mm；

和/或，根系吸水模型子单元，用于构建并调整根系吸水模型：

根系吸水量的表达式为：

$S=T_p\&times;\frac{e^{-AROOT\&times;z_2}-e^{-AROOT\&times;z_1}}{e^{-AROOT\&times;LR}-1},$

其中T_p为植被潜在蒸腾量，AROOT为描述根主要分布的参数，z₁、z₂分别为所求土壤层垂直方向上的两端坐标，单位m，LR为根系深度，单位m；

和/或，层间土壤水分增量模型子单元，用于构建并调整层间土壤水分增量模型：

当进入第i层的剩余入渗量大于AD_max,i时，降水或灌溉导致的层间土壤水分增量：AD_i＝AD_max,i，其中AD_max,i为各层的最大增量，AD_max,i＝(W_f,i-W_i)·L_i，i为土壤层数，L_i为土层厚度，取10mm；

否则，层间土壤水分增量等于剩余入渗水量，公式为：

$Gd=\left\{\begin{array}{c}F-\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i(F>\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i)\\ 0\end{array}\right.,$

其中，F＝P+I-S_V-Rs，P为期内降水量，单位mm，I为灌水量，在非灌区及天然植被区为0，单位mm；S_V为期内累计植被截留量，单位mm；Rs为期内地表径流量，单位mm；

和/或，根系层补给量模型子单元，用于构建并调整根系层补给量模型：

因上层根系吸水损失的水分表示公式为：

S_uppersop_loss_i＝K_si*ET_p*(number_layer-i)/depth_sum

其中，number_layer＝10，depth_sum＝100,单位为mm。

进一步地，所述蒸散发模型子单元包括：

参考作物蒸散量计算模块，用于计算参考作物蒸散量ET_p0，公式为：

${\mathrm{ET}}_{p0}=\frac{0.408\mathrm{\&Delta;}(R_n-G)+\mathrm{\&gamma;}\frac{900}{T_a+273}{U}_2(e_S-e_d)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+0.34U_2)},$

式中：△是气温T时的饱和水汽压曲线斜率，单位kPa·℃^-1：

T是日平均温度，单位℃；

γ是干湿表常数，单位kPa·℃^-1：C_p是空气定压比热，取1.013×10^-3MJ·kg^-1·℃^-1；P_r是大气压，单位kPa；λ是蒸发潜热，取2.45MJ·kg^-1；ε是水汽分子量与干空气分子量之比，取0.622；

e_S是气温T下的饱和水汽压，单位kPa：

e_d是实际水汽压，单位kPa：e_d＝RH×e_s，RH是空气相对湿度，单位％；

U₂是2米高度处风速，单位m·s^-1：U_z为z高度处观测到的风速，单位m·s^-1；

R_n是净辐射，单位MJ·m^-2·d^-1：

$R_n=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_n\left(t\right)dt=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_{n\_max}sin\&lsqb;\left(\frac{t-t_{rise}}{t_{set}-t_{rise}}\right)\mathrm{\&pi;}\&rsqb;dt,$

t_rise，t_set分别为日出时间和日落时间，R_{n_max}为日最大净辐射，t为卫星过境时刻的当地时间；

G是土壤热通量，单位MJ·m^-2·d^-1：有植被覆盖的地面为：

G＝R_n[Γ_c+(1-f_c)(Γ_s-Γ_c)]

式中：全植被覆盖下，土壤热通量与净辐射的比值Γ_c＝0.05；裸地情况下，土壤热通量与净辐射比值Γ_s＝0.315，f_c为植被覆盖率；

对于水体和冰雪，G＝0.5R_n；

作物系数计算模块，用于计算作物系数K_c：

非生长期的K_c小于0.4；

生长期的K_c为：K_c＝K_cb+0.05，K_cb为非完全覆盖条件下的基本植物系数：

式中，K_cbmin为在缺少植被情况下的裸地最小植物系数，取0.15-0.20；K_cbfull为全覆盖条件下，叶面积指数大于3的基本植物系数；h为植被高度，单位m；f_cell为有效植被覆盖度。

进一步地，

所述模型参数包括：气象参数、下垫面植被参数、下垫面土壤参数、经验参数和其它参数中的一种或几种；

所述气象参数包括：降雨量、灌溉量、气温、地表温度、2米处风速、大气压、空气相对湿度、生长期内平均风速和平均最小相对湿度中的一种或几种；

所述下垫面植被参数包括：叶面积指数、植被盖度、地表反照率、波段31发射率、波段32发射率、植被高度和平均植被最大高度中的一种或几种；

所述下垫面土壤参数包括：初始土壤含水量、土壤田间含水量、土壤饱和含水量、土壤萎蔫含水量中的一种或几种；

所述经验参数包括：地表径流参数1、地表径流参数2和地表覆被因子中的一种或几种；

所述其它参数包括：中分辨率成像光谱仪过境时间、日出时间、日落时间和地理纬度中的一种或几种。

(三)有益效果

可见，在本发明提出的一种多层土壤水分模拟方法和系统中，能够将遥感反演方法和传统的典型土壤水分层均衡模型相耦合，构建出物理机制明确、结构简单的模型，并能保持一定的土壤水模拟精度，实现多层土壤水分模拟的空间化和动态化。

本发明充分考虑了区域尺度的地表不均一性,将空间离散化为像元，在像元尺度上获取传统的多层土壤水模型中各子模块所需的部分参数，建立耦合遥感信息的多层土壤水模拟模型与方法，进行降水植被截留、土壤含水量、径流、蒸散发、层间水分增量和根系补给量的全过程模拟，实现了对多层土壤含水量的实时监控，从而有利于加强区域水资源管理，合理利用降水资源，缓解水资源匮乏。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例多层土壤水分模拟方法的基本流程示意图；

图2是本发明一个优选实施例多层土壤水分模拟方法的模拟流程图；

图3是本发明一个优选实施例多层土壤水分模拟方法在渭河流域的模拟值与实测值之间的相关关系图；

图4为是本发明一个优选实施例多层土壤水分模拟方法在2008年渭河流域进行土壤水分模拟的空间分布图；

图5是本发明实施例多层土壤水分模拟系统的基本结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例首先提供一种多层土壤水分模拟方法，参见图1，包括：

步骤101：构建土壤水分层均衡模型，结合遥感技术调整模型结构。

步骤102：利用遥感技术获取模型参数，构建流域水文空间信息数据库。

步骤103：利用调整后的土壤水分层均衡模型和所述流域水文空间信息数据库，开展多层土壤水分过程数值模拟。

可见，在本发明实施例提出的一种多层土壤水分模拟方法中，能够将遥感反演方法和传统的典型土壤水分层均衡模型相耦合，构建出物理机制明确、结构简单的模型，并能保持一定的土壤水模拟精度，实现多层土壤水分模拟的空间化和动态化。

土壤水分层均衡模型相对简单，在一般情况下也能够达到一定的模拟精度，是模拟和预报土壤水分变化最有效的方法之一。本发明实施例考虑地下水位较深的情况，忽略地下水对根系层的补给作用，考虑植被截留和地表径流过程；在水平方向上，只对有植被覆盖区进行模拟。在本发明的一个实施例中，优选地，土壤水分层均衡模型可以包括：土壤水分平衡模型、植被截流模型、地表径流模型、蒸散发模型、根系吸水模型、层间土壤水分增量模型、根系层补给量模型中的一个或多个。

随着遥感技术的发展，空间变异性问题得到了有效解决，本发明实施例将遥感技术引入到土壤水分层均衡模型数值模拟中。土壤水分层均衡模型的遥感驱动可从模型结构、参数获取方法两个方面出发，以土壤水分层均衡模型为基础框架，在继承模型优点的同时，增强部分模型参数的物理意义，同时增加模型对遥感数据的耦合性。在实际应用中，可以采用易于获取的遥感数据与产品为模型的主要输入，构建遥感驱动的土壤水分层均衡模型，以期达到削弱土壤水分模拟对地面观测数据的依赖性的目的。在利用遥感驱动开展土壤水分层均衡模型中部分子模块结构的改进，加强模型对遥感数据的耦合的过程中，可以在蒸散发模型中参考作物蒸散发选用Penman-Monteith公式进行获取；还可以考虑植被截留和地表径流，其中植被截留模型选择Hoyningen-Huene，该模型直接建立截留量与植被特征参数(叶面积指数、植被盖度)之间的关系，模型参数均能通过遥感技术获取，地表径流模型采用地表径流采用VGTM模型计算；也可以改进根系吸水模型，增强与遥感技术的耦合性。

优选地，可以具体进行如下操作：

1、构建并调整土壤水分平衡模型：土壤水分平衡是指一定时间内，作物根系层范围一定深度内水分的收支关系。水量交换侧重垂直方向，根据渭河流域地下水位情况，在研究中忽略地下水对根系层的补给作用。若考虑时段初根系层土壤贮水量为W₁，时段末为W₂，根系层土壤贮水量的变化用下式表示：

W₂-W₁＝(P+I)-(S_V+R_S+E_a+T_a+D)

式中，P为期内降水量，单位mm；I为灌水量，在非灌区及天然植被区为0，单位mm；S_V为期内累计植被截留量，单位mm；Rs为期内地表径流量，单位mm；E_a为期内累计土壤蒸发量，单位mm；T_a为期内累计植被蒸腾量，单位mm；D为根系层底部水分向下渗漏项，单位mm。

若将土壤根区分为若干层次，则(P+I)-(S_V+R)转变为各层水分增量ΔD，植被蒸腾量转变为各层根系吸水量S，无降水条件下，各层水分净上移量表示为Δs_upperrop_loss，则分层土壤水分平衡方程表示为：

降水条件下：

(表层，有棵间蒸发)

$W_i^{j+1}=W_i^j+{\mathrm{\&Delta;D}}_i^j-S_i^j(i\&GreaterEqual;2)$

非降水条件下：

(表层)

(i≥2)

Δs_upperrop_loss_i＝s_upperrop_loss_i+1-s_upperrop_loss_i-1(i≥2)

式中，i为土壤层次，j为时间节点，E_a为棵间蒸发量(发生在表层)，ΔD_i为降水或灌水后第i层土壤水分增量，S_i为第i层作物根系的吸水量，s_upperrop_loss_i为第i层净水分上移量。

2、构建并调整植被截留模型：

植被截留是指大气降水到达冠层后，部分降水被植被的冠层(树干和枝叶)截留并存储的现象，它对雨水具有在数量和时间上重新分配的功能，截留水量将以蒸发的形式返回大气中，并影响同期的蒸散发能力。我国学者对地跨我国南北不同气候带及相应的森林植被类型林冠截留率的分析研究表明，截留率变动范围在11.4％～34.3％，变动系数6.68％～55.5％，可见林冠截留对水文过程的影响。林冠截留模型包括经验模型、半经验半理论模型以及理论模型，其中以Rutter模型和Gash解析模型较为完善且被广泛应用，但是模型的参数较难获取，因此，本发明实施例采用半经验半理论模型。

降雨过程中植被对降水的截留主要有枝叶作用，因此截留量主要与植被特征有关，如植被盖度、叶面积、叶面粗糙度，随植被类型和生长时段而变化。Aston(1979)分析林冠截留和植被盖度之间的关系，并构建了林冠截留量计算方程：

$S_v=c_v\&times;S_{max}\&times;\&lsqb;1-e^{-\mathrm{\&eta;}\frac{P}{S_{max}}}\&rsqb;$

式中：S_v为累计截留量，单位mm；c_v为植被盖度，单位％，反映覆被空间分布情况；P为累积降水量，单位mm；S_max为树冠蓄水能力，单位mm，即林冠最大截留量；η为校正系数。林冠最大截留量主要取决于叶面特征(叶面积、叶面粗糙度)。Hoyningen-Huene(1981)建立了直接基于LAI估算最大截留量的公式：

S_max＝0.935+0.498×LAI-0.00575×LAI²

η＝0.046×LAI

在以上植被截留模型中，模型输入参数包括降雨量、叶面积指数(LAI)和植被盖度，遥感技术在这些参数的获取上均已有较为成熟的应用，分辨率成像光谱仪(moderate-resolution imaging spectroradiometer，MODIS)等提供较高时间分辨率的LAI产品，能为模型提供驱动数据。

3、构建并调整地表径流模型：

降雨经林冠截留后到达表层土壤，一部分在表面发生地表径流，剩余部分进入表层土壤补充土壤含水量。地表径流采用VGTM模型计算，模型认为地表径流与有效降雨呈非线性关系，这种关系通过时变增益因子表示，时变增益因子通过表层土壤水分求算，同时增加下垫面覆被状况对产流状况的影响通过地表覆被因子C来表达，其主要与植被生长状况有关。地表径流模型计算公式如下：

$Rs=g_1{\left(\frac{{\mathrm{AW}}_u}{{\mathrm{WM}}_u\&CenterDot;C}\right)}^{g_2}\&CenterDot;P^{\&prime;}$

式中：Rs为地表产流量，单位mm；AW_u为表层土壤湿度，单位％；WM_u为表层土壤饱和含水量，单位％；P'为有效降雨量，单位mm，为降雨量扣除林冠截留后到达地面的净雨量；g₁与g₂是时变增益因子的有关参数(0<g₁<1，1<g₂)，其中g₁为土壤饱和后径流系数，g₂为土壤水影响系数；C为覆被影响参数。

4、构建并调整蒸散发模型：

实际蒸散发模型：

实际蒸散量(ET_a)取决于潜在蒸散量(ET_p)以及制约蒸散过程的土壤水分胁迫系数(K_S)。因此，实际蒸散发可表示为：

ET_a＝K_s×ET_p

Jensen等(1972)的研究成果表明，土壤水分胁迫系数K_s用下式表示：

K_s＝ln(A_v+1)/ln(101)

A_v＝[(W-W_m)/(W_f-W_m)]×100

式中，W_f为田间持水量，单位mm；W_m为萎蔫系数，单位mm；A_v为相对有效含水量。

因为植被覆盖区实际蒸散量由土表实际棵间蒸发量(E_a)和实际蒸腾量(T_a)组成，则上式可改写成：

E_a＝K_s×E_p

T_a＝K_s×T_p

对于土壤棵间蒸发，水分胁迫系数K_s计算过程中，W_m＝W_a，W_a为风干土含水量。土表潜在棵间蒸发量(E_p)由Ritchie(1972)公式求得：

$E_p=\left\{\begin{array}{cc}{\mathrm{ET}}_p\&times;(1-0.43\&times;LAI)& LAI\&le;1\\ {\mathrm{ET}}_p/1.1\&times;\exp (-0.4\&times;LAI)& LAI>1\end{array}\right.$

植被潜在蒸腾量(T_p)表示为

T_p＝ET_p-E_p

潜在蒸散发模型：

潜在蒸散发是指蒸发面在特定的气象条件下充分供水时的蒸发量或蒸发率，对某一蒸发面来说，其蒸发能力就是在同样气象条件下可能达到的最大蒸发量，气象站和水文站一般采用D20和E601蒸发皿观测。国内外也提出了大量潜在蒸散发估算方法，大致可分为气温法、辐射法、综合法、水量平衡法和质量守恒法等。

本文采用FAO推荐的方法计算潜在蒸散(Allen,1998)，潜在蒸散的表达式如下

ET_p＝K_c·ET_p0

式中：K_c为作物系数，通过遥感估算获得；ET_p0是参考作物蒸散量(mm)，通过遥感估算获得。

参考作物蒸散ET_p0：

计算参考作物蒸散量ET_p0的方法很多，归纳起来大致可分为经验公式法、水汽扩散法、能量平衡法、综合法和BP神经网络模型等几大类。为使计算公式统一化、标准化，联合国粮农组织FAO在1998年出版的《Crop evapotranspiration，Guidelines for computingcrop water requirements》一书中推荐用修正的Penman-Monteith公式作为计算参考作物蒸发蒸腾量ET_p0的标准化方法。

${\mathrm{ET}}_{p0}=\frac{0.408\mathrm{\&Delta;}(R_n-G)+\mathrm{\&gamma;}\frac{900}{T_a+273}{U}_2(e_S-e_d)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+0.34U_2)}$

式中：△是气温T时的饱和水汽压曲线斜率(kPa·℃^-1)；R_n是净辐射(MJ·m^-2·d^-1)；G是土壤热通量(MJ·m^-2·d^-1)；γ是干湿表常数(kPa·℃^-1)；U₂是2米处风速(m·s^-1)；e_S是气温T下的饱和水汽压(kPa)；e_d是实际水汽压(kPa)。

A.△饱和水汽压曲线斜率：

$\mathrm{\&Delta;}=\frac{4098\&lsqb;0.6108\exp \left(\frac{17.27T}{T+237.3}\right)\&rsqb;}{{(T+237.3)}^2}$

式中：△是气温T时的饱和水汽压曲线斜率(kPa·℃^-1)；T是日平均温度(℃)。

B.干湿表常数：

$\mathrm{\&gamma;}=\frac{C_p{P}_r}{\mathrm{\&epsiv;}\mathrm{\&lambda;}}=0.665\&times;{10}^{-3}{P}_r$

式中：γ是干湿表常数(kPa·℃^-1)；C_p是空气定压比热，值为1.013×10^-3MJ·kg^-1·℃^-1(指一定气压下，单位体积的空气温度升高1℃所需的能量为1.013×10^-3MJ·kg^-1·℃^-1)；P_r是大气压(kPa)；λ是蒸发潜热，取2.45MJ·kg^-1；ε是水汽分子量与干空气分子量之比，值为0.622。

C.平均饱和水汽压e_S：

$e_s=0.6108\exp \left(\frac{17.27T}{T+237.3}\right);$

D.实际水汽压e_d：

e_d＝RH×e_s

式中：RH是空气相对湿度(％)。

E.风速：

$U_2=\frac{4.87}{ln(67.8z-5.42)}{U}_z$

式中：U_z为z高度处观测到的风速(m·s^-1)；U₂为2m高度处的风速。

F.地表净辐射：

地表净辐射通量(Net Radiation)又称辐射平衡(Radiation Balance)或辐射差额(Radiation Budget)，是指地表面净得的短波辐射与长波辐射的和，即地表辐射能量收支的差额。它是地表面能量、动量、水分输送与交换过程中的主要能源，是地气能量交换中最重要的组分，也是驱动地表蒸散和显热通量的主要来源，近年来，基于MODIS估算晴空条件下地表净辐射及分量(短波净辐射、长波净辐射)的研究较多(Gautam，2005；Youngryel，2008；Tang，2008；Wang，2009)。本文基于能量平衡法，采用MODIS数据估算卫星过境时刻的瞬时地表净辐射：

INR＝R_S ^↓-R_S ^↑+R_L ^↓-R_L ^↑＝(1-α)R_S+R_L ^↓-R_L ^↑

式中：R_S ^↓，R_S ^↑，R_L ^↓，R_L ^↑分别为下行、上行短波辐射和下行、上行长波辐射，w/m²；α为地表反照率；R_S为向下的短波辐射，即到达地表的太阳总辐射。

长波辐射采用Steffan–Boltzmann方程：

R_L ^↓-R_L ^↑＝σε_aε_sT_a ⁴-σε_sT_s ⁴

式中：ε_a为大气发射率，ε_s为地表发射率，T_a为空气温度，单位K；T_s为地表温度，单位K；σ为Steffan–Boltzmann常数，取值5.67×10^-8W/m².K^-4。

空气比辐射率ε_a计算公式如下：

ε_a＝9.2×10^-6×T_a ²

地表发射率ε_s采用Liang(2004)的非线性公式计算：

ε_s＝0.273+1.778ε₃₁-1.807ε₃₁ε₃₂-1.037ε₃₂+1.774ε² ₃₂

式中：ε₃₁，ε₃₂分别为波段31和波段32的发射率。

采用类似于Lagouarde and Brunet(1983)用正弦函数求算地表温度的方法，以求算净辐射的日过程：

$R_n\left(t\right)=R_{n\_\max }sin\&lsqb;\left(\frac{t-t_{rise}}{t_{set}-t_{rise}}\right)\mathrm{\&pi;}\&rsqb;$

式中：R_{n_max}为日最大净辐射，t为卫星过境时刻的当地时间，t_rise，t_set分别对应为净辐射变为正和变为负的当地时间，分别为日出时间和日落时间。

日最大净辐射可由瞬时净辐射求得：

$R_{n\_max}=\frac{INR}{sin\&lsqb;\left(\frac{t_{pass}-t_{rise}}{t_{set}-t_{rise}}\right)\&rsqb;}$

由此，日净辐射总量计算公式为：

$R_n=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_n\left(t\right)dt=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_{n\_max}sin\&lsqb;\left(\frac{t-t_{rise}}{t_{set}-t_{rise}}\right)\mathrm{\&pi;}\&rsqb;dt,$

日均净辐射计算公式如下：

$DANR=\frac{2INR}{\mathrm{\&pi;}sin\&lsqb;\left(\frac{T-2a}{2T}\right)\mathrm{\&pi;}\&rsqb;}$

式中：T为日长，a为净辐射达到最大R_{n_max}的时间和MODIS过境时间的时差。R_{n_max}一般出现在当地时间的12：30，MODIS过境时间Time由遥感产品获得。

日长为日落时间t_set与日出时间t_rise的时差，由下式求算：

T＝12+[2arcsin(tgL·tgδ)]/15°

式中：L为纬度，δ为太阳赤纬。

G.土壤热通量：

净辐射能量一部分被用作蒸散发，一部分被用来加热地表的大气，剩下的则被储藏在土壤或水体之中，这部分土壤或水体的热交换能量即土壤热通量。目前土壤热通量还不能直接通过遥感技术估算，一般通过它与净辐射及植被覆盖或者植被指数的关系来确定，本文采用Su(2002)提出的计算方法，对于有植被覆盖的地面，计算公式如下：

G＝R_n[Γ_c+(1-f_c)(Γ_s-Γ_c)]

式中：全植被覆盖下，土壤热通量与净辐射的比值Γ_c＝0.05；裸地情况下，土壤热通量与净辐射比值Γ_s＝0.315，f_c为植被覆盖率。

对于水体和冰雪，土壤热通量的计算多采用与R_n取比值的方法，本文水体、积雪与冰川的土壤热通量与净辐射的比值取0.5(Waters.R.，et al，2002)，即认为进入水体和冰雪的净辐射有一半被吸收成为冰雪的G：G＝0.5R_n。

植物系数K_c：

自然植被生长不如农田那么密集，在稀疏条件下(叶面积指数LAI<3)时，植物系数K_c受到叶面积指数(或者覆盖度)的影响，FAO专门提出了适合于自然植被非完全覆盖条件下植物系数的计算方法。

K_c包括生长期和非生长期两种。非生长期的K_c一般都小于0.4，如果有地表覆盖，参数取值则要高一些。生长期的K_c可采用下式计算：

K_c＝K_cb+0.05

式中，K_c为生长期植物系数，K_cb为非完全覆盖条件下的基本植物系数，可由下式计算：

$K_{cb}=K_{cbmin}+(K_{cbfull}-K_{cbmin})\&times;min(1,2\&times;f_{cell}^{1/(1+h)})$

式中，K_cbmin为在缺少植被情况下的裸地最小植物系数(0.15-0.20)；K_cbfull为全覆盖条件下(LAI>3)基本植物系数；h为植被高度(m)；f_cell为有效植被覆盖度，可采用下式计算：

f_cell＝f_c/sinη

式中，f为植被覆盖度，取0.01-1；sinη为太阳高于海平面的中间角正弦值，在最大蒸散法期间一般在11.0-15.0。f_cell通常取中午12点的计算值，此时sinη的计算公式为：

式中，δ和分别为太阳赤纬(rad)、纬度(rad)。

$K_{cbfull}=K_{cb,h}+\&lsqb;0.04(U_2-2)-0.004({\mathrm{RH}}_{min}-45)\&rsqb;\&times;{\left(\frac{h}{3}\right)}^{0.3}$

式中，K_cb，h为在标准湿润和风速条件下(RH_min＝45％和U₂＝2m/s)的全覆盖植被的K_cb基本植物系数；U₂为生长期内2m高位置处的平均风速；RH_min为生长期的平均最小相对湿度(％)；h为平均植被最大高度(m)。K_cb,h可以利用下式进行估计：

K_cb,h＝1.0+0.1×h h≤2m

K_cb,h＝1.2 h>2m

5、构建并调整根系吸水模型：

有关作物的根系吸水函数的研究很多，但许多模型的形式复杂，需要很多参数。有很多研究使用DeJong吸水函数，在本研究中在DeJong吸水函数结构基础上，相对根密度采用Kristensen-Jensen模型中根系分布函数一项，对每一层的根系吸水量进行估算。

原模型中根系吸水量表达式为：

$S_{mi}=\frac{T_p*{\mathrm{RD}}_i}{{\&Integral;}_1^N{\mathrm{RD}}_{i}dz}$

S_i＝K_si*S_mi

式中：S_mi为第i层内作物根系吸收的水分；RD_i为根密度。

Kristensen-Jensen模型中根系分布函数表达式为：

$RDF=\frac{{\&Integral;}_{z_1}^{z_2}R\left(z\right)dz}{{\&Integral;}_0^{LR}R\left(z\right)dz}$

且：

logR(z)＝logR₀-AROOT×z

RDF的计算形式可转化为与植被根系深度有关的形式：

$RDF=\frac{e^{-AROOT\&times;z_2}-e^{-AROOT\&times;z_1}}{e^{-AROOT\&times;LR}-1}$

式中，RDF为根系分布函数，AROOT为描述根主要分布的参数，z₁，z₂分别为所求土壤层垂直方向上的两端坐标，单位m；LR为根系深度，单位m。

本研究中根系吸水量表达式转化为如下形式：

$S=T_p\&times;\frac{e^{-AROOT\&times;z_2}-e^{-AROOT\&times;z_1}}{e^{-AROOT\&times;LR}-1}.$

6、构建并调整层间土壤水分增量模型：

简化土壤入渗水再分配过程，不考虑分配时间，分配为依次使各层达到田间持水量后多余水量下渗，假定各层的最大增量为AD_max,i，则有：

AD_max,i＝(W_f,i-W_i)·L_i

式中：i为土壤层数；L_i为土层厚度，在本文中为10mm。当进入i层的剩余入渗量大于AD_max,i时，

降水或灌溉导致的土壤水分增量AD_i＝AD_max,i，否则AD_i等于剩余入渗水量。

土壤水渗漏按下式计算：

$Gd=\left\{\begin{array}{c}F-\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i(F>\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i)\\ 0\end{array}\right.,$

式中，F＝P+I-S_V-Rs。

7、构建并调整根系层补给量模型：

在实际情况下，随着主要根系层内根系吸水，土层之间的水势梯度逐渐增加，下层土壤内的水分不断地向上补给，这部分水的损失与蒸散有关，据经验，由于上层根系吸水损失的水分表示为：

S_uppersop_loss_i＝K_si*ET_p*(number_layer-i)/depth_sum

在本发明实施例中，number_layer＝10，depth_sum＝100,单位为mm。在本发明的一个实施例中，遥感驱动可以充分利用公用平台提供的遥感和GIS数据产品，为模型提供下垫面参数，这些参数均可从公用数据平台直接或者间接地获取，从而构建流域水文空间信息数据库，进而完成模型与遥感数据的耦合，开展流域土壤水分过程数值模拟。本发明实施例中，在渭河流域构建的土壤水分均衡模型的各个子模块中涉及的参数及其获取方法见表1，模型的参数大致可分为4类：气象参数、下垫面植被参数、下垫面土壤参数和各类经验参数。其中大部分气象参数和下垫面植被参数可基于公用数据平台的遥感和GIS数据产品确定，下垫面土壤参数结合相应的土壤数据库进行确定，模型涉及的部分经验参数可结合已有相关研究成果确定，有些参数则有待于率定。

表1模型参数及其确定方法

本发明实施例中，多层土壤水分模拟流程见图2。在本发明的一个实施例中，对渭河流域多层土壤水分模拟的模拟值与实测值之间的相关关系图见图3，其中，R2代表监测值和模拟值的相关系数的平方值，代表两者的相关性大小；土壤水分的单位为mm。可见，本发明实施例方法的监测值和模拟值之间具有很好的对应关系。图4为利用本发明实施例方法，在2008年渭河流域进行土壤水分模拟的空间分布图。

本发明一个实施例还提供了一种多层土壤水分模拟系统，参见图5，包括：

模型构建单元501，用于构建土壤水分层均衡模型，结合遥感技术调整模型结构；

遥感参数单元502，用于利用遥感技术获取模型参数，构建流域水文空间信息数据库；

模拟单元503，用于利用调整后的土壤水分层均衡模型和所述流域水文空间信息数据库，开展多层土壤水分过程数值模拟。

在本发明的一个实施例中，优选地，模型构建单元501可以包括：

模型分类子单元，用于将土壤水分层均衡模型分为：土壤水分平衡模型、植被截流模型、地表径流模型、蒸散发模型、根系吸水模型、层间土壤水分增量模型、根系层补给量模型中的一种或多种类型。

在本发明的另一个实施例中，优选地，模型构建单元501可以包括：

土壤水分平衡模型子单元，用于构建并调整土壤水分平衡模型：

降水条件下，表层、有棵间蒸发的土壤水分平衡方程为：

$W_1^{j+1}=W_1^j+{\mathrm{\&Delta;D}}_1^j-S_1^j,$

其余为：

非降水条件下，表层土壤水分平衡方程为：

$W_1^{j+1}=W_1^j+s\_uppersop\_{\mathrm{loss}}_2-{\mathrm{Ea}}^j-S_1^j,$

其余为：

Δs_upperrop_loss_i＝s_upperrop_loss_i+1-s_upperrop_loss_i-1(i≥2)；

其中i为土壤层次，j为时间节点，为第i层时段初根系层土壤贮水量，为第i层时段末根系层土壤贮水量，E_a为表层棵间蒸发量，ΔD_i为降水或灌水后第i层土壤水分增量，S_i为第i层作物根系的吸水量，s_upperrop_loss_i为第i层净水分上移量；

和/或，植被截流模型子单元，用于构建并调整植被截流模型：

构建林冠截留量计算方程为：

$S_v=c_v\&times;S_{max}\&times;\&lsqb;1-e^{-\mathrm{\&eta;}\frac{P}{S_{max}}}\&rsqb;$

式中，S_v为累计截留量，单位mm；c_v为植被盖度，单位％，反映覆被空间分布情况；P为累积降水量，单位mm；S_max为树冠蓄水能力，单位mm，即林冠最大截留量；η为校正系数；

其中林冠最大截留量为：

S_max＝0.935+0.498×LAI-0.00575×LAI²，

校正系数：η＝0.046×LAI，LAI为叶面积指数；

和/或，地表径流模型子单元，用于构建并调整地表径流模型：

$Rs=g_1{\left(\frac{{\mathrm{AW}}_u}{{\mathrm{WM}}_u\&CenterDot;C}\right)}^{g_2}\&CenterDot;P^{\&prime;}$

式中，Rs为地表产流量，单位mm；AW_u为表层土壤湿度，单位％；WM_u为表层土壤饱和含水量，单位％；P'为有效降雨量，单位mm，是降雨量扣除林冠截留后到达地面的净雨量；g₁与g₂是时变增益因子的有关参数，0<g₁<1，1<g₂，其中g₁为土壤饱和后径流系数，g₂为土壤水影响系数；C为覆被影响参数；

和/或，蒸散发模型子单元，用于构建并调整蒸散发模型：

植被覆盖区实际蒸散量由土表实际棵间蒸发量E_a和实际蒸腾量T_a组成，公式为：

E_a＝K_s×E_p，

T_a＝K_s×T_p，

其中土壤水分胁迫系数K_s为：

K_s＝ln(A_v+1)/ln(101)，

A_v＝[(W-W_a)/(W_f-W_a)]×100，

式中，W_f为田间持水量，单位mm；W_a为风干土含水量；

土表潜在棵间蒸发量E_p为：

LAI为叶面积指数，

植被潜在蒸腾量T_p为：

T_p＝ET_p-E_p；

潜在蒸散ET_p的表达式为：

ET_p＝K_c·ET_p0，

式中，K_c为作物系数；ET_p0为参考作物蒸散量，单位mm；

和/或，根系吸水模型子单元，用于构建并调整根系吸水模型：

根系吸水量的表达式为：

$S=T_p\&times;\frac{e^{-AROOT\&times;z_2}-e^{-AROOT\&times;z_1}}{e^{-AROOT\&times;LR}-1},$

其中T_p为植被潜在蒸腾量，AROOT为描述根主要分布的参数，z₁、z₂分别为所求土壤层垂直方向上的两端坐标，单位m，LR为根系深度，单位m；

和/或，层间土壤水分增量模型子单元，用于构建并调整层间土壤水分增量模型：

当进入第i层的剩余入渗量大于AD_max,i时，降水或灌溉导致的层间土壤水分增量：AD_i＝AD_max,i，其中AD_max,i为各层的最大增量，AD_max,i＝(W_f,i-W_i)·L_i，i为土壤层数，L_i为土层厚度，取10mm；

否则，层间土壤水分增量等于剩余入渗水量，公式为：

$Gd=\left\{\begin{array}{c}F-\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i(F>\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{ADm}}_i)\\ 0\end{array}\right.,$

其中，F＝P+I-S_V-Rs，P为期内降水量，单位mm，I为灌水量，在非灌区及天然植被区为0，单位mm；S_V为期内累计植被截留量，单位mm；Rs为期内地表径流量，单位mm；

和/或，根系层补给量模型子单元，用于构建并调整根系层补给量模型：

因上层根系吸水损失的水分表示公式为：

S_uppersop_loss_i＝K_si*ET_p*(number_layer-i)/depth_sum

其中，number_layer＝10，depth_sum＝100,单位为mm。

在本发明的一个实施例中，优选地，蒸散发模型子单元可以包括：

参考作物蒸散量计算模块，用于计算参考作物蒸散量ET_p0，公式为：

${\mathrm{ET}}_{p0}=\frac{0.408\mathrm{\&Delta;}(R_n-G)+\mathrm{\&gamma;}\frac{900}{T_a+273}{U}_2(e_S-e_d)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+0.34U_2)},$

式中：△是气温T时的饱和水汽压曲线斜率，单位kPa·℃^-1：

T是日平均温度，单位℃；

γ是干湿表常数，单位kPa·℃^-1：C_p是空气定压比热，取1.013×10^-3MJ·kg^-1·℃^-1；P_r是大气压，单位kPa；λ是蒸发潜热，取2.45MJ·kg^-1；ε是水汽分子量与干空气分子量之比，取0.622；

e_S是气温T下的饱和水汽压，单位kPa：

e_d是实际水汽压，单位kPa：e_d＝RH×e_s，RH是空气相对湿度，单位％；

U₂是2米高度处风速，单位m·s^-1：U_z为z高度处观测到的风速，单位m·s^-1；

R_n是净辐射，单位MJ·m^-2·d^-1：

$R_n=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_n\left(t\right)dt=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_{n\_max}sin\&lsqb;\left(\frac{t-t_{rise}}{t_{set}-t_{rise}}\right)\mathrm{\&pi;}\&rsqb;dt,$

t_rise，t_set分别为日出时间和日落时间，R_{n_max}为日最大净辐射，t为卫星过境时刻的当地时间；

G是土壤热通量，单位MJ·m^-2·d^-1：有植被覆盖的地面为：

G＝R_n[Γ_c+(1-f_c)(Γ_s-Γ_c)]

式中：全植被覆盖下，土壤热通量与净辐射的比值Γ_c＝0.05；裸地情况下，土壤热通量与净辐射比值Γ_s＝0.315，f_c为植被覆盖率；

对于水体和冰雪，G＝0.5R_n；

作物系数计算模块，用于计算作物系数K_c：

非生长期的K_c小于0.4；

生长期的K_c为：K_c＝K_cb+0.05，K_cb为非完全覆盖条件下的基本植物系数：

$K_{cb}=K_{cbmin}+(K_{cbfull}-K_{cb\min })\&times;\min (1,2\&times;f_{cell}^{1/(1+h)}),$

式中，K_cbmin为在缺少植被情况下的裸地最小植物系数，取0.15-0.20；K_cbfull为全覆盖条件下，叶面积指数大于3的基本植物系数；h为植被高度，单位m；f_cell为有效植被覆盖度。

在本发明的一个实施例中，优选地，模型参数可以包括：气象参数、下垫面植被参数、下垫面土壤参数、经验参数和其它参数等，其中：

气象参数可以包括：降雨量、灌溉量、气温、地表温度、2米处风速、大气压、空气相对湿度、生长期内平均风速和平均最小相对湿度等。

下垫面植被参数可以包括：叶面积指数、植被盖度、地表反照率、波段31发射率、波段32发射率、植被高度和平均植被最大高度等。

下垫面土壤参数可以包括：初始土壤含水量、土壤田间含水量、土壤饱和含水量、土壤萎蔫含水量等。

经验参数可以包括：地表径流参数1、地表径流参数2和地表覆被因子等。

其它参数可以包括：中分辨率成像光谱仪过境时间、日出时间、日落时间和地理纬度等。

可见，本发明实施例具有如下有益效果：

在本发明实施例提出的一种多层土壤水分模拟方法和系统中，能够将遥感反演方法和传统的典型土壤水分层均衡模型相耦合，构建出物理机制明确、结构简单的模型，并能保持一定的土壤水模拟精度，实现多层土壤水分模拟的空间化和动态化。

本发明实施例充分考虑了区域尺度的地表不均一性,将空间离散化为像元，在像元尺度上获取传统的多层土壤水模型中各子模块所需的部分参数，建立耦合遥感信息的多层土壤水模拟模型与方法，进行降水植被截留、土壤含水量、径流、蒸散发、层间水分增量和根系补给量的全过程模拟，实现了对多层土壤含水量的实时监控，从而有利于加强区域水资源管理，合理利用降水资源，缓解水资源匮乏。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (6)

1.一种多层土壤水分模拟方法，其特征在于，包括：

构建包括蒸散发模型的土壤水分层均衡模型，结合遥感技术调整土壤水分层均衡模型的结构，包括：构建并调整蒸散发模型：

植被覆盖区实际蒸散量由表层棵间蒸发量E_a和实际蒸腾量T_a组成：

E_a＝K_s×E_p，

T_a＝K_s×T_p，

其中土壤水分胁迫系数K_s为：

K_s＝ln(A_v+1)/ln(101)，

A_v＝[(W-W_a)/(W_f-W_a)]×100，

式中，W为根系层土壤贮水量，W_f为田间持水量，单位mm；W_a为风干土含水量，单位mm；A_v为相对有效含水量；

土表潜在棵间蒸发量E_p为：

LAI为叶面积指数，

植被潜在蒸腾量T_p为：

T_p＝ET_p-E_p；

潜在蒸散ET_p的表达式为：

ET_p＝K_c·ET_p0，

式中，K_c为作物系数；ET_p0为参考作物蒸散量，单位mm；

其中，所述参考作物蒸散量ET_p0的计算方法为：

${\mathrm{ET}}_{p0}=\frac{0.408\mathrm{\&Delta;}(R_n-G)+\mathrm{\&gamma;}\frac{900}{T_a+273}{U}_2(e_S-e_d)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+0.34U_2)},$

式中：△是气温T时的饱和水汽压曲线斜率，单位kPa·℃^-1

T是日平均温度，单位℃；

γ是干湿表常数，单位kPa·℃^-1：C_p是空气定压比热，取1.013×10^-3MJ·kg^-1·℃^-1；P_r是大气压，单位kPa；λ是蒸发潜热，取2.45MJ·kg^-1；ε是水汽分子量与干空气分子量之比，取0.622；

e_S是气温T下的饱和水汽压，单位kPa：

e_d是实际水汽压，单位kPa：e_d＝RH×e_S，RH是空气相对湿度，单位％；

U₂是2米高度处风速，单位m·s^-1：U_z为z高度处观测到的风速，单位m·s^-1；

R_n是净辐射，单位MJ·m^-2·d^-1：

$R_n=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_n\left(t\right)dt=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_{n\_max}sin\&lsqb;\left(\frac{t-t_{rise}}{t_{set}-t_{rise}}\right)\mathrm{\&pi;}\&rsqb;dt,$

t_rise，t_set分别为日出时间和日落时间，R_{n_max}为日最大净辐射，t为卫星过境时刻的当地时间；

G是土壤热通量，单位MJ·m^-2·d^-1：有植被覆盖的地面为：

G＝R_n[Γ_c+(1-f_c)(Γ_s-Γ_c)]

式中：全植被覆盖下，土壤热通量与净辐射的比值Γ_c＝0.05；裸地情况下，土壤热通量与净辐射比值Γ_s＝0.315，f_c为植被覆盖率；

对于水体和冰雪，G＝0.5R_n；

所述作物系数K_c的计算方法为：

非生长期的K_c小于0.4；

生长期的K_c为：K_c＝K_cb+0.05，K_cb为非完全覆盖条件下的基本植物系数：

$K_{cb}=K_{cb\min }+(K_{cbfull}-K_{cb\min })\&times;min(1,2\&times;f_{cell}^{1/(1+h)}),$

式中，K_{cb min}为在缺少植被情况下的裸地最小植物系数，取0.15-0.20；K_cbfull为全覆盖条件下，叶面积指数大于3的基本植物系数；h为植被高度，单位m；f_cell为有效植被覆盖度；

利用遥感技术获取土壤水分层均衡模型参数，构建流域水文空间信息数据库；

利用调整后的土壤水分层均衡模型和所述流域水文空间信息数据库，开展多层土壤水分过程数值模拟。

2.根据权利要求1所述的多层土壤水分模拟方法，其特征在于，所述土壤水分层均衡模型还包括：土壤水分平衡模型、植被截流模型、地表径流模型、根系吸水模型、层间土壤水分增量模型、根系层补给量模型中的一个或多个。

3.根据权利要求1或2所述的多层土壤水分模拟方法，其特征在于：

所述土壤水分层均衡模型参数包括：气象参数、下垫面植被参数、下垫面土壤参数、经验参数和其它参数；

所述气象参数包括：大气压、2米处风速、空气相对湿度、气温，还包括：降雨量、灌溉量、地表温度、生长期内平均风速和平均最小相对湿度中的一种或几种；

所述下垫面植被参数包括：叶面积指数、植被高度，还包括：植被盖度、地表反照率、波段31发射率、波段32发射率和平均植被最大高度中的一种或几种；

所述下垫面土壤参数包括：土壤田间含水量，还包括：初始土壤含水量、土壤饱和含水量、土壤萎蔫含水量中的一种或几种；

所述经验参数包括：地表径流参数1、地表径流参数2和地表覆被因子中的一种或几种；

所述其它参数包括：日出时间、日落时间，还包括：中分辨率成像光谱仪过境时间和地理纬度中的一种或几种。

4.一种多层土壤水分模拟系统，其特征在于，包括：模型构建单元、遥感参数单元和模拟单元，模型构建单元和遥感参数单元分别与模拟单元相连，其中：

模型构建单元，用于构建包括蒸散发模型的土壤水分层均衡模型，结合遥感技术调整土壤水分层均衡模型的结构，

包括：蒸散发模型子单元，用于构建并调整蒸散发模型：

植被覆盖区实际蒸散量由表层棵间蒸发量E_a和实际蒸腾量T_a组成，公式为：

E_a＝K_s×E_p，

T_a＝K_s×T_p，

其中土壤水分胁迫系数K_s为：

K_s＝ln(A_v+1)/ln(101)，

A_v＝[(W-W_a)/(W_f-W_a)]×100，

式中，W为根系层土壤贮水量，W_f为田间持水量，单位mm；W_a为风干土含水量，单位mm；A_v为相对有效含水量；

土表潜在棵间蒸发量E_p为：

LAI为叶面积指数，

植被潜在蒸腾量T_p为：

T_p＝ET_p-E_p；

潜在蒸散ET_p的表达式为：

ET_p＝K_c·ET_p0，

式中，K_c为作物系数；ET_p0为参考作物蒸散量，单位mm；

其中，所述蒸散发模型子单元包括：

参考作物蒸散量计算模块，用于计算参考作物蒸散量ET_p0，公式为：

${\mathrm{ET}}_{p0}=\frac{0.408\mathrm{\&Delta;}(R_n-G)+\mathrm{\&gamma;}\frac{900}{T_a+273}{U}_2(e_S-e_d)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+0.34U_2)},$

式中：△是气温T时的饱和水汽压曲线斜率，单位kPa·℃^-1：

T是日平均温度，单位℃；

γ是干湿表常数，单位kPa·℃^-1：C_p是空气定压比热，取1.013×10^-3MJ·kg^-1·℃^-1；P_r是大气压，单位kPa；λ是蒸发潜热，取2.45MJ·kg^-1；ε是水汽分子量与干空气分子量之比，取0.622；

e_S是气温T下的饱和水汽压，单位kPa：

e_d是实际水汽压，单位kPa：e_d＝RH×e_S，RH是空气相对湿度，单位％；

U₂是2米高度处风速，单位m·s^-1：U_z为z高度处观测到的风速，单位m·s^-1；

R_n是净辐射，单位MJ·m^-2·d^-1：

$R_n=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_n\left(t\right)dt=\underset{t_{rise}}{\overset{t_{set}}{\&Integral;}}{R}_{n\_max}sin\&lsqb;\left(\frac{t-t_{rise}}{t_{set}-t_{rise}}\right)\mathrm{\&pi;}\&rsqb;dt,$

t_rise，t_set分别为日出时间和日落时间，R_{n_max}为日最大净辐射，t为卫星过境时刻的当地时间；

G是土壤热通量，单位MJ·m^-2·d^-1：有植被覆盖的地面为：

G＝R_n[Γ_c+(1-f_c)(Γ_s-Γ_c)]

式中：全植被覆盖下，土壤热通量与净辐射的比值Γ_c＝0.05；裸地情况下，土壤热通量与净辐射比值Γ_s＝0.315，f_c为植被覆盖率；

对于水体和冰雪，G＝0.5R_n；

作物系数计算模块，用于计算作物系数K_c：

非生长期的K_c小于0.4；

生长期的K_c为：K_c＝K_cb+0.05，K_cb为非完全覆盖条件下的基本植物系数：

$K_{cb}=K_{cb\min }+(K_{cbfull}-K_{cb\min })\&times;min(1,2\&times;f_{cell}^{1/(1+h)}),$

式中，K_{cb min}为在缺少植被情况下的裸地最小植物系数，取0.15-0.20；K_cbfull为全覆盖条件下，叶面积指数大于3的基本植物系数；h为植被高度，单位m；f_cell为有效植被覆盖度；

遥感参数单元，用于利用遥感技术获取土壤水分层均衡模型参数，构建流域水文空间信息数据库；

模拟单元，用于利用调整后的土壤水分层均衡模型和所述流域水文空间信息数据库，开展多层土壤水分过程数值模拟。

5.根据权利要求4所述的多层土壤水分模拟系统，其特征在于，所述模型构建单元包括：

模型分类子单元，用于将土壤水分层均衡模型分为：蒸散发模型、土壤水分平衡模型、植被截流模型、地表径流模型、根系吸水模型、层间土壤水分增量模型、根系层补给量模型。

6.根据权利要求4或5所述的多层土壤水分模拟系统，其特征在于：

所述土壤水分层均衡模型参数包括：气象参数、下垫面植被参数、下垫面土壤参数、经验参数和其它参数；

所述气象参数包括：大气压、2米处风速、空气相对湿度、气温，还包括：降雨量、灌溉量、地表温度、生长期内平均风速和平均最小相对湿度中的一种或几种；

所述下垫面植被参数包括：叶面积指数、植被高度，还包括：植被盖度、地表反照率、波段31发射率、波段32发射率和平均植被最大高度中的一种或几种；

所述下垫面土壤参数包括：土壤田间含水量，还包括：初始土壤含水量、土壤饱和含水量、土壤萎蔫含水量中的一种或几种；

所述经验参数包括：地表径流参数1、地表径流参数2和地表覆被因子中的一种或几种；

所述其它参数包括：日出时间、日落时间，还包括：中分辨率成像光谱仪过境时间和地理纬度中的一种或几种。