一种农业区水资源环境诊断方法及其系统.pdf

本发明涉及一种农业区水资源环境诊断方法，其特征在于，所述方法包括：S1，建立蒸散发模型和产流模型，根据所述蒸散发模型和所述产流模型对农业区用水进行诊断；S2，建立非点源污染模型，对农业区非点源污染进行诊断；S3，根据所述农业区用水的诊断和所述非点源污染的诊断，对农业区水资源环境进行诊断。本发明的方法既考虑不同种植结构的下垫面非均一性带来的总的用水和灌溉用水空间差异，又考虑了不同下垫面非点源污染排放强度，能够有效地诊断水资源环境问题。

1.  一种农业区水资源环境诊断方法，其特征在于，所述方法包括：
S1，建立蒸散发模型和产流模型，根据所述蒸散发模型和所述产流模型对农业区用水进行诊断；
S2，建立非点源污染模型，对农业区非点源污染进行诊断；
S3，根据所述农业区用水的诊断和所述非点源污染的诊断，对农业区水资源环境进行诊断。

2.  根据权利要求1所述的一种农业区水资源环境诊断方法，其特征在于，所述步骤S1具体包括：
S11，建立蒸散发模型，计算农业区的蒸散发量；
S12，建立产流模型，计算农业区的径流量；
S13，根据所述蒸散发量和所述径流量，计算农业区用水诊断值S_t；
S14，若S_t＜1,则说明农业区用水量在控制的范围内，反之则超用。

3.  根据权利要求2所述的一种农业区水资源环境诊断方法，其特征在于，在所述步骤S11具体包括：
S111，计算植被截留蒸发量：
E_can＝min(S_v,ET_P',P)
其中，S_v为植被截留量，Sv=cv×Smax×[1-e-ηPcumSmax],]]>
c_v为植被覆盖度，P_cum为累积降水量，S_max为林冠最大截留量，且S_max＝0.935+0.498×LAI-0.00575×LAI²；η为校正系数，η＝0.046×LAI，ET_P'为潜在蒸散发量，P为降水量，LAI为叶面积指数；
S112，计算植被蒸腾量：
E_at＝f₁(LAI)×f₂(θ)×(RDF)×(ET_P'-E_can)
其中，RDF为根系分布函数；f₁(LAI)为蒸腾对植被叶面积的依赖度，f₂(θ)为根系层土壤水分函数；(ET_P'-E_can)为扣除林冠截留量后的最大蒸发能力。
S113，计算土壤蒸发量：
E_s＝ET_P'×f₃(θ)+(ET_P'-E_at-ET_P'×f₃(θ))×f₄(θ)×(1-f₁(LAI))
其中，f₃(θ)为，f₄(θ)为；
S114，计算所述农业区的蒸散发量：
ET＝E_can+E_at+E_s。

4.  根据权利要求2所述的一种农业区水资源环境诊断方法，其特征在于，所述步骤S12具体包括：
S121，计算地表径流量：
Rs=g1(AWuWMu·C)g2·P′,]]>
其中，AW_u为表层土壤湿度，WM_u为表层土壤饱和含水量，P'为有效降雨量，g₁为土壤饱和后径流系数，g₂为土壤水影响系数0<g₁<1，1<g₂；C为覆被影响参数；
S122，计算表层壤中流量：
Rss＝AW_u·K_r·Thick，
其中，AW_u为表层土壤含水量，K_r为土壤水出流系数；Thick为表层土壤厚度；
S123，计算深层壤中流：
Rds＝AW_d·K_d·Thick_s
其中，AW_d为深层土壤含水量，K_d为深层土壤径流出流系数，Thick_s为深层土壤厚度；
S124，计算地下径流：
Rg＝AW_g·K_g·Thick_g
其中，AW_g为地下层含水量，K_g为地下径流出流系数，Thick_g为地下层土壤厚度；
S125，计算所述农业区的径流量：R＝Rs+Rss+Rsd+Rg。

5.  根据权利要求2所述的一种农业区水资源环境诊断方法，其特征在于，所述步骤S13具体包括：
S131，建立降水平衡公式：
P＝ET_p+R+ΔAW，
式中，ET_p为降水中用于蒸散发部分的量，R为径流量，ΔAW为土壤湿度增加量；
S132，计算农业区总用水量：
W＝W_I+W_p
其中，W_I＝ET-ET_P＝ET+R-P，W_I为灌溉用水量，ET为农业区的蒸散量，W_p为因种植污染等级为劣Ⅲ类水量，
S133，根据所述农业区总用水量计算农业种植用水诊断值：
S_t＝(W_I+W_p)/W_s
其中，W_s为可供用水量。

6.  根据权利要求1所述的一种农业区水资源环境诊断方法，其特征在于，所述步骤S2具体包括：
S21，计算溶解态非点源污染的值：
C=Σi=14&epsiv;&epsiv;0&CenterDot;(1-ekR0t)&CenterDot;Qi&CenterDot;Ni&CenterDot;Si]]>
其中，C为污染负荷，i为溶解态污染类型，n为类型数；Q_i为污染源强；k为地面冲刷系数；R₀为标准雨强,ε为径流系数，ε₀为标准径流系数，N_i为自然因子修正系数；S_i为社会因子修正系数；
S22，计算吸附态非点源污染的值：
C_a＝X·Q_a·τ
式中,X为土壤侵蚀量；Q_a为流失土壤中的氮磷污染物背景含量；τ为吸附态污染物富集系数；
S23，计算非点源污染量的值：
NP＝C+C_a
S24，计算子区域水环境的诊断值：S_l＝NP_l/R_l，其中，l为自然数，表示第l个子区域；
S25，计算所述农业区的水环境的诊断值：
S＝∑S_l。

7.  根据权利要求6所述的一种农业区水资源环境诊断方法，其特征在于，所述污染源强的计算如下：
S 71，计算农业生产的污染源强：
Qa=ΣxnΣymAx&CenterDot;ωxy&CenterDot;M]]>
其中，x为农田类型，n为农田类型数，y为农田污染物指标类型，m为农田染物指标类型的类型数，A为农田面积，ω_xy为各种作为不同污染物指标的源强系数，M为修正系数；
S72，计算畜禽养殖污染源强：
Qb=ΣhnnΣkmmCh&CenterDot;Dh&CenterDot;kh&CenterDot;ωhk&CenterDot;ηhk,]]>
其中，h为畜禽养殖类型，nn为畜禽养殖类型数，k为畜禽养殖污染物指标类型，mm为畜禽养殖污染物指标类型数，C为畜禽数量，D为不同畜禽类型饲养周期，k_h为不同畜禽类型粪便排泄指数，ω_hk为不同畜禽类型不同污染物指标的源强系数；η_hk为不同畜禽类型不同污染物指标的流失率；
S73，计算农村居民点的污染源强：
Qc=ΣfmnPc&CenterDot;ωf,]]>
其中，f为污染物指标类型，mn为类型数，P_c为乡村人口数量，ω_f为不同污染物指标的源强系数；
S74，计算城镇居民点的污染源强：
Qd=ΣgnmPd&CenterDot;ωg,]]>
其中，g为污染物指标类型，nm为类型数，P_d为城镇人口数量，ω_g为不同污染物指标的源强系数；
S75，计算所述农业区的污染源强：Q＝Qa+Qb+Qc+Qd。

8.  一种农业区水资源环境诊断系统，其特征在于，上述系统包括：
农业区用水诊断模块，用于建立蒸散发模型和产流模型，根据所述蒸散发模型和所述产流模型对农业区用水进行诊断；
农业区非点源污染诊断模块，用于建立非点源污染模型，对农业区非点源污染进行诊断；
农业区水资源环境诊断模块，根据所述农业区用水的诊断和所述非点源污染的诊断，对农业区水资源环境进行诊断。

一种农业区水资源环境诊断方法及其系统
技术领域
本发明涉及生态水文和遥感技术领域，具体地，涉及一种农业区水资源环境的诊断方法及其系统。
背景技术
农业区指以种植业为主，既生产粮食作物与经济作物，同时兼营林、牧、渔各业的地区。是粮食和经济作物产品的最大提供者，也是轻工业原料、畜禽产品、水产品等各种农副产品的主要产地。
我国干旱和半干旱地区面积占到了国土面积的近60％，而农业用水仍延续粗放式管理，据统计，我国农田单位用水的粮食产量不足2.4斤/立方米，而世界先进水平为4斤/立方米左右。随着近些年地表与地下水资源持续超采导致生态环境恶化、地下水位下降等问题十分突出。与此同时，农业区农田、居民点、畜禽养殖区等农药、化肥过量施用、废弃与排泄物无序排放产生的大量非点源污染物严重污染着水环境，直接威胁者农村饮用水、灌溉和生态环境用水安全，实行严格的水环境管理制度，推进农业灌溉用水总量控制和定额管理势在必行。
农业区的种植用水主要来自大气降水和灌溉用水，其中，大气降水到地面后分流成冠层截留、补充土壤水和径流三个部分，而土壤水又以土壤蒸发和植被蒸腾(ET)的形式耗散到大气中；灌溉用水主要用于补充土壤水，同样以土壤蒸发和植被蒸腾的形式耗散到大气中。
传统的地表蒸散估算方法局限于单点尺度，而且大范围密集观测成本较高。对于区域而言，下垫面情况复杂且非均一，土地类型和小气候环境差异大，传统地表蒸散单点监测和计算方法局限性日 益突出。遥感技术的出现和发展为这个问题的解决带来了新的希望，可见光、近红外和热红外等波段能够提供与水分和能量平衡过程密切相关的一些参数。蒸散的测定可分为直接测量与间接估算，直接测定地表蒸散发主要采用的是蒸发皿和蒸发池、蒸渗仪、波文比-能量平衡观测系统、涡动相关仪(EC)、大孔径闪烁仪(LAS)等。间接测定或估算方法有水量平衡法、Penman-Montieth公式、Priestley-Taylor模型、互补相关方法、作物缺水指数法等。目前，利用遥感监测地表蒸散量的方法，概括起来可分为：(1)经验统计方法；(2)与传统计算方法相结合模型；(3)基于能量平衡方程模型；(4)数值模拟以及数据同化方法。
非点源污染是相对点源污染而言，是指溶解的和固体的污染物从非特定的地点，在降水(或融雪)冲刷作用下，通过径流过程而汇入受纳水体(包括河流、湖泊、水库和海湾等)并引起水体的富营养化或其它形式的污染(Novotny和Olem，1993)。美国清洁水法修正案(1997)对非点源污染的定义为：污染物以广域的、分散的、微量的形式进入地表及地下水体。
非点源污染根据污染物产生和迁移转化过程可分为溶解态非点源污染和吸附态非点源污染。溶解态非点源污染具有水溶性，污染物伴随水文产汇流过程聚集；吸附态非点源污染是指污染物吸附在土壤颗粒上，污染物汇集多受土壤侵蚀过程控制。基于上述两种污染物产生机理，将非点源污染负荷分为溶解态非点源污染负荷和吸附态非点源污染负荷分别进行计算。相关非点源污染模型包括SWAT、AWIMM、新安江等国内外模型，近些年我国郝芳华等(2006)针对我国的大尺度区域特征，结合自然因子和社会因子建立了二元结构溶解态非点源污染模型。该模型综合了物理模型和统计模型的优势，比较适合我国地形地貌特征的非点源污染负荷估算。杨胜天等(2006)以USLE土壤侵蚀方程为基础，建立了吸附态氮和磷流 失量的估算模型，并进行了全国境内水土流失影响下吸附态氮磷流失量估算。
上述的蒸散发模型(ET)和非点源污染估算模型分别以能量或水平衡过程和以径流过程为非点源污染负荷产生的基本动力过程进行蒸散发和非点源污染估算的，结合遥感和GIS空间技术可以对下垫面进行农业种植用水和农业区所产生的非点源污染进行估算。然而，系统的针对农业区种植结构进行种植用水，尤其是灌溉用水估算存在明显不足，很难揭示超采地下或过量使用河流湖泊水量。同时，对于非点源污染除则更多的是污染负荷量的估算，因为缺乏径流量分析而难以得出单位径流量中所含的污染物量(即评估单元非污染排放强度＝非点源污染负荷/径流量)的诊断依据。
故存在，针对不同农业区的种植结构，不能有系统效地解析出农业区不同种植结构的下垫面灌溉用水量和评估单元非污染排放强度的空间格局，面向农业区精细化水资源环境管理具有较大的局限性。
发明内容
本发明针对上述情况，提供一种农业区水资源环境诊断方法及系统，该方法既考虑不同种植结构的下垫面非均一性带来的总的用水和灌溉用水差异，又考虑了不同下垫面非点源污染排放强度，能够有效地诊断水资源环境问题。
根据本发明的一个方面，提供一种农业区水资源环境诊断方法，所述方法包括：S1，建立蒸散发模型和产流模型，根据所述蒸散发模型和所述产流模型对农业区用水进行诊断；S2，建立非点源污染模型，对农业区非点源污染进行诊断；S3，根据所述农业区用水的诊断和所述非点源污染的诊断，对农业区水资源环境进行诊断。
其中，所述步骤S1具体包括：S11，建立蒸散发模型，计算农业区的蒸散发量；S12，建立产流模型，计算农业区的径流量；S13，根 据所述蒸散发量和所述径流量，计算农业区用水诊断值S_t；S14，若S_t＜1,则说明农业区用水量在控制的范围内，反之则超用。
其中，在所述步骤S11具体包括：S111，计算植被截留蒸发量：
E_can＝min(S_v,ET_P',P)
其中，S_v为植被截留量，Sv=cv×Smax×[1-e-ηPcumSmax],]]>
c_v为植被覆盖度，P_cum为累积降水量(单位：mm)，S_max为林冠最大截留量(单位：mm)，且S_max＝0.935+0.498×LAI-0.00575×LAI²；η为校正系数，η＝0.046×LAI，ET_P'为潜在蒸散发量(单位：mm)，P为降水量(单位：mm)，LAI为叶面积指数；
S112，计算植被蒸腾量：E_at＝f₁(LAI)×f₂(θ)×(RDF)×(ET_P'-E_can)
其中，RDF为根系分布函数；f₁(LAI)为蒸腾对植被叶面积的依赖度，f₂(θ)为根系层土壤水分函数；(ET_P'-E_can)为扣除林冠截留量后的最大蒸发能力。
S113，计算土壤蒸发量：
E_s＝ET_P'×f₃(θ)+(ET_P'-E_at-ET_P'×f₃(θ))×f₄(θ)×(1-f₁(LAI))
其中，f₃(θ)为，f₄(θ)为；
S114，计算所述农业区的蒸散发量：ET＝E_can+E_at+E_s
其中，所述步骤S12具体包括：S121，计算地表径流量：
Rs=g1(AWuWMu&CenterDot;C)g2&CenterDot;P′,]]>
其中，AW_u为表层土壤湿度(单位：mm)，WM_u为表层土壤饱和含水量(单位：mm)，P'为有效降雨量(单位：mm)，g₁为土壤饱和后径流系数，g₂为土壤水影响系数0<g₁<1，1<g₂；C为覆被影响参数；
S122，计算表层壤中流量：Rss＝AW_u·K_r·Thick，其中，AW_u为表层土壤含水量(单位：mm)，K_r为土壤水出流系数；Thick为表层土壤厚度；
S123，计算深层壤中流：Rds＝AW_d·K_d·Thick_s，其中，AW_d为深层土壤含水量(单位：mm)，K_d为深层土壤径流出流系数，Thick_s为深层土壤厚度；
S124，计算地下径流：Rg＝AW_g·K_g·Thick_g，其中，AW_g为地下层含水量(单位：mm)，K_g为地下径流出流系数，Thick_g为地下层土壤厚度；
S125，计算所述农业区的径流量：R＝Rs+Rss+Rsd+Rg。
其中，所述步骤S13具体包括：S131，建立降水平衡公式：
P＝ET_p+R+ΔAW，式中，ET_p为降水中用于蒸散发部分的量(单位：mm)，R为径流量(单位：mm)，ΔAW为土壤湿度增加量(单位：mm)；S132，计算农业区总用水量：
W＝W_I+W_p，其中，W_I＝ET-ET_P＝ET+R-P，W_I为灌溉用水量(单位：mm)，ET为农业区的蒸散量(单位：mm)，W_p为因种植污染等级为劣Ⅲ类水量，
S133，根据所述农业区总用水量计算农业种植用水诊断值：
S_t＝(W_I+W_p)/W_s，其中，W_s为可供用水量(单位：mm)。
其中，所述步骤S2具体包括：S21，计算溶解态非点源污染的值：
C=Σi=14&epsiv;&epsiv;0&CenterDot;(1-ekR0t)&CenterDot;Qi&CenterDot;Ni&CenterDot;Si]]>
C为污染负荷，i为溶解态污染类型，n为类型数；Q_i为污染源强；k为地面冲刷系数；R₀为标准雨强,ε为径流系数，ε₀为标准径流系数，N_i为自然因子修正系数；S_i为社会因子修正系数；
S22，计算吸附态非点源污染的值：C_a＝X·Q_a·τ，式中,式中,X为土壤侵蚀量；Q_a为流失土壤中的氮磷污染物背景含量；τ为吸附态污染物富集系数；
S23，计算非点源污染量的值：NP＝C+C_a
S24，计算子区域水环境的诊断值：S_l＝NP_l/R_l，其中，l为自然数，表示第l个子区域；
S25，计算所述农业区的水环境的诊断值：S＝∑S_l。
其中，所述污染源强的计算如下：
S71，计算农业生产的污染源强：
其中，x为农田类型，n为农田类型数，y为农田污染物指标类型，m为农田染物指标类型的类型数，A为农田面积，ω_xy为各种作为不同污染物指标的源强系数，M为修正系数；
S72，计算畜禽养殖污染源强：
其中，h为畜禽养殖类型，nn为畜禽养殖类型数，k为畜禽养殖污染物指标类型，mm为畜禽养殖污染物指标类型数，C为畜禽数量，D为不同畜禽类型饲养周期，k_h为不同畜禽类型粪便排泄指数，ω_hk为不同畜禽类型不同污染物指标的源强系数；η_hk为不同畜禽类型不同污染物指标的流失率；
S73，计算农村居民点的污染源强：
其中，f为污染物指标类型，mn为类型数，P_c为乡村人口数量，ω_f为不同污染物指标的源强系数；
S74，计算城镇居民点的污染源强：其中，g为污染物指标类型，nm为类型数，P_d为城镇人口数量，ω_g为不同污染物指标的源强系数；
S75，计算所述农业区的污染源强：Q＝Qa+Qb+Qc+Qd。
根据本发明的另一个方面，提供一种农业区水资源环境诊断系统，上述系统包括：
农业区用水诊断模块，用于建立蒸散发模型和产流模型，根据所述蒸散发模型和所述产流模型对农业区用水进行诊断；
农业区非点源污染诊断模块，用于建立非点源污染模型，对农业区非点源污染进行诊断；
农业区水资源环境诊断模块，根据所述农业区用水的诊断和所述非点源污染的诊断，对农业区水资源环境进行诊断。
通过上述实施例可知，本发明的一种农业区水资源环境的诊断方法及其系统，基于蒸散发模型、径流模型和非点源污染模型，以年为时间尺度，并综合了径流数据，有效地解决了农业区不同种植结构总用水和灌溉用水估算问题和非点源污染输出强度估算问题。同时解决了传统非遥感技术估算方法中因未考虑下垫面差异的问题，并最终能够有效地诊断水资源环境问题。
附图说明
通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，在附图中：
图1示出了本发明的一种农业区水资源环境诊断的方法的流程图。
图2示出了本发明的一种农业区水资源环境诊断系统的结构框图。
具体实施方式
下面将结合附图对本发明的实施例进行详细描述。
图1示出了本发明的一种农业区水资源环境诊断的方法的流程图。
参照图1，本发明实施例的一种农业区水资源环境诊断方法包括步骤：
S1、建立蒸散发模型和产流模型，根据所述蒸散发模型和所述产流模型对农业区用水进行诊断。
具体包括步骤：
S11，建立蒸散发模型，计算农业区的蒸散发量；
在本实施例中，蒸散发部分的计算使用以下方法。
蒸散发包括植被截留蒸发、土壤蒸发、水面蒸发和植被蒸腾 等。对于水域的蒸发比较简单，可认为其等同于潜在蒸散发量；而对于非水面的蒸散法，需考虑植被截留蒸发、土壤蒸发和植被蒸腾3个方面。
S111，计算植被截留蒸发量
植被截留是指大气降水到达冠层后，部分降水被植被的冠层(树干和枝叶)截留并存储的现象，它对雨水具有在数量和时间上重新分配的功能，截留水量将以蒸发的形式返回大气中，并影响同期的蒸散发能力。我国学者对地跨我国南北不同气候带及相应的森林植被类型林冠截留率的分析研究表明，截留率变动范围在11.4％～34.3％，变动系数6.68％～55.5％，可见林冠截留对水文过程的影响。
林冠截留模型包括经验模型、半经验半理论模型以及理论模型，其中以Rutter模型和Gash解析模型较为完善且被广泛应用，但是模型的参数较难获取，基于此，本实施例采用半经验半理论模型。
降雨过程中植被对降水的截留主要由枝叶作用，因此截留量主要与植被特征有关，如植被盖度、叶面积、叶面粗糙度，随植被类型和生长时段而变化。因此经过分析林冠截留和植被盖度之间的关系，构建了林冠截留量计算方程：
Sv=cv×Smax×[1-e-ηPcumSmax]]]>
式中：S_v为植被截留量；c_v为植被覆盖度，反映覆被空间分布情况；P_cum为累积降水量；S_max为树冠蓄水能力即林冠最大截留量；η为校正系数；林冠最大截留量主要取决于叶面特征(叶面积、叶面粗糙度)。
最大截留量的公式为：
S_max＝0.935+0.498×LAI-0.00575×LAI²
η＝0.046×LAI
其中，LAI为叶面积指数。
在以上植被截留模型中，模型输入参数包括降雨量、叶面积指数和植被盖度，遥感技术在这些参数的获取上均已有较为成熟的应用，MODIS、VEGETATION等提供较高时间分辨率的LAI产品，能为模型提供驱动数据。
模型将冠层截留过程模拟为截留贮水过程，并认为植被截留的水分最终均以蒸发的形式返回到大气中。因此，蒸散发优先计算植被截留蒸发，其大小则取决于植被截留量和潜在蒸散发能力，取两者中的低值，若冠层截留有足够的水分，则林冠截留蒸发等于潜在蒸散发大小，公式如下：
E_can＝min(S_v,ET_P',P)
式中：E_can为植被截留蒸发量(单位：mm)；ET_P'为潜在蒸散发量(单位：mm)，可通过遥感数据估算；P为降水量(单位：mm)。
S112，计算植被蒸腾量
植物蒸腾速率主要取决于土壤的水分状况和植被覆被状况。对于植被状况，模型中引入LAI和根系深度表达植被对蒸腾作用的影响，并将植被蒸腾表达为LAI、根系区土壤水分以及根系密度的函数。当土壤水分充足达到田间持水量时，植被的蒸腾作用不受土壤水分含量的影响，仅与植被特征相关；而当土壤含水量达到萎蔫含水量，即植物最大水分输出量小于可能蒸腾量时，植物气孔自行关闭，蒸腾作用停止。模型表达式如下：
E_at＝f₁(LAI)×f₂(θ)×(RDF)×(ET_P'-E_can)
式中：E_at为植被蒸腾量(单位：mm)；RDF为根系分布函数；f₁和RDF反映植被覆盖对植被蒸腾的影响，分别表示蒸腾对植被叶面积和根系密度的依赖度；f₂为根系层土壤水分函数，反映土壤水分状况对蒸腾的影响；(ET_P'-E_can)为扣除林冠截留量后的最大蒸发能力。f₁、f₂和RDF的表达式如下：
f₁(LAI)＝max(0,min(1,(C₁+C₂·LAI)))
f2(θ)=0θ≤θw1-(θf-θθf-θw)C3ETpθw≤θ≤θf1θ&GreaterEqual;θf]]>
且：RDF=&Integral;z1z2R(z)dz&Integral;0LRR(z)dz]]>
log R(z)＝log R₀-AROOT×z
RDF的计算形式可转化为与植被根系深度有关的形式：
RDF=e-AROOT×z2-e-AROOT×z1e-AROOT×LR-1]]>
式中，θ为土壤含水量(单位：mm)；θ_f为田间持水量(单位：mm)；θ_w为萎蔫含水量(单位：mm)；C₁，C₂，C₃为参数；AROOT为描述根主要分布的参数，z1，z2分别为所求土壤层垂直方向上的两端坐标；L_R为根系深度。
S113，计算土壤蒸发量
Kristensen-Jensen模型认为土壤蒸发仅发生在表层土壤，由非饱和上层土壤基础蒸发量和土壤达到田间持水量时多余的水分蒸发组成，其中土壤蒸发受表层土壤水分的限制，同时取决于潜在蒸散发，当土壤含水量降至剩余土壤含水量时，土壤水蒸发停止。土壤蒸发计算公式如下：
Es＝ET_p·f₃(θ)+(ET_p-E_at-ET_p·f₃(θ))·f₄(θ)·(1-f₁(LAI))
f3(θ)=0θ≤θrC2(θ/θw)θr<θ<θwC2θ&GreaterEqual;θw]]>
f4(θ)=θ-0.5(θw+θf)θf-0.5(θw+θf)θ&GreaterEqual;0.5(θw+θf)0θ<0.5(θw+θf)]]>
S114，计算农业区的总的蒸散发量。
非水面单元上的蒸散发为植被截留蒸发、植被蒸腾、土壤蒸发之和。模型的输入参数主要包括决定蒸发能力的潜在蒸散发、影响植被截留和植被蒸腾的LAI、影响植被蒸腾的根系深度、影响植被蒸腾和土壤蒸发的土壤水分特征参数。
基于上述，本实施例的真散发部分为上述植被截留蒸发、植被蒸腾、土壤蒸发之和，即：
ET＝E_can+E_at+E_s
S12，建立产流模型，计算农业区的径流量；
具体步骤包括：
产流模型是DTVGM的核心模块，其通过时变增益因子优先计算地表径流，并耦合水量平衡方程，以牛顿迭代的方法求算土壤湿度，再求算各个水文分量。目前，遥感在地表径流和壤中流计算中的应用较少，因此，对于地表及地表以下的产流过程，RS-DTVGM和原模型保持一致，地表径流采用TVGM，壤中流采用自由水蓄水库线性出流。
S121，计算地表径流量。
降雨经林冠截留后到达表层土壤，一部分在表面发生地表径流，剩余部分进入表层土壤补充土壤含水量。地表径流采用VGTM模型计算，模型认为地表径流与有效降雨呈非线性关系，这种关系通过时变增益因子表示，时变增益因子通过表层土壤水分求算，同时增加下垫面覆被状况对产流状况的影响通过地表覆被因子C来表达，其主要与植被生长状况有关。地表径流模型计算公式如下：
Rs=g1(AWuWMu&CenterDot;C)g2&CenterDot;P′]]>
式中：AW_u为表层土壤湿度(单位：mm)；WM_u为表层土壤饱和含水量(单位：mm)；P'为有效降雨量(单位：mm)，为降雨 量扣除林冠截留后到达地面的净雨量和融雪量之和；g₁与g₂是时变增益因子的有关参数(0<g₁<1，1<g₂)，其中g₁为土壤饱和后径流系数，g₂为土壤水影响系数；C为覆被影响参数。
S122，计算表层壤中流
扣除地表径流外的降雨以一定的速率下渗进入土壤，一部分保存在土壤中，使土壤含水量增大；一部分沿着土壤空隙流动，从坡侧土壤空隙流出，转换为地表径流，注入河槽的部分径流称为表层壤中流，采用自由水蓄水库线性出流计算，公式如下：
Rss＝AW_u·K_r·Thick
式中：AW_u为表层土壤含水量(单位：mm)；K_r为土壤水出流系数；Thick为表层土壤厚度。在实际计算时土壤湿度取时段起止土壤湿度的平均值：
AWu=AWui+AWui+12]]>
式中：AWu_i和AWu_i+1分别为时段开始和结束时的土壤水含量。
S123，计算深层壤中流与地下径流
表层土壤水分在重力和水势作用下向深层下渗，给定表层到深层的下渗率为fc，即可求算上层土壤渗入到下层的水量，地下径流亦相同。
深层壤中流和地下径流计算方式同于表层壤中流，即
Rds＝AW_d·K_d·Thick_s
Rg＝AW_g·K_g·Thick_g
式中：Rds和Rg分别为深层壤中流和地下径流(单位：mm)；AW_d和AW_g分别为深层土壤和地下层含水量(单位：mm)；K_d和K_g分别深层土壤和地下径流出流系数，地下径流一般出流较为稳定，在冰期，径流的主要供给为地下流，因此可根据冬季径流量确定出流系数。
S125，计算所述农业区的总的径流量：
总的径流量为地表径流、壤中流、地下径流之和：
R＝Rs+Rss+Rsd+Rg
S13，根据蒸散发量和径流量，计算农业区用水诊断值。
根据水平衡原理，农业区用水包括农业种植、居民点、养殖等用水，其中居民生活用水、养殖用水等可以结合统计分析等方式获得用水量，而农业区种植用水则需要考虑降水、径流、灌溉用水、蒸散发等信息，所以是用水模型的核心。
对农业种植用水估算如下过程。
建立降水平衡公式：
P＝ET_p+R+ΔAW
式中，P为降水量(单位：mm)，是降雨强度和持续时长的积；ET_p为降水中用于蒸散发的部分(单位：mm)，R为径流量(单位：mm)，ΔAW为土壤湿度增加量(单位：mm)，当以年为估算单位时，这里ΔAW忽略不计，取近似值0。所以，ET_p可以近似为：
ET_p＝P-R
对于总得蒸散发量的计算公式为：ET＝P-R+W_I
式中，ET为种植区(可以是子区域)总的蒸散发量(单位：mm)，W_I为灌溉用水量(单位：mm)，进一步，灌溉用水量可以表达为；
W_I＝ET-ET_P＝ET+R-P
则总的用水量的计算公式为：W＝W_I+W_p式中，W_p为因种植污染等级为劣Ⅲ类水量。
因此，根据农业区总用水量计算农业种植用水诊断值：
S_t＝(W_I+W_p)/W_s
W_s可供用水量(单位：mm)。
S14，诊断是否用水过量，若S_t＜1,则说明农业区用水量在控制的范围内，反之则超用。
在上述实施例中，蒸散发模型和产流模型均需要土壤含水量作为输入，虽然目前基于遥感反演土壤含水量已有相当多的研究，但仍不能直接应用于水文模型中作为模型驱动。光学反演土壤水分多是根据地表覆被建立土壤水分与遥感参数(如NDVI等)之间的统计关系，多用于较长时间段的土壤水分含量的估算；微波方法探测土壤含水量方面，被动微波空间分辨率太粗，探测的为表层很薄的一层土壤，对于深层的无法获得，而主动微波成本高，在区域尺度应用还不具可行性。因此，选择结合水量平衡方程，采用牛顿跌代法求算土壤含水量。
对于表层土壤，水量平衡方程如下：
P_i+AWu_i＝AWu_i+1+Rs_i+ETa_i+Rss_i+WUB_i
式中：P_i为时段内降雨量(单位：mm)；AWu_i为时段初表层土壤含水量(单位：mm)；AWu_i+1为时段末表层土壤含水量(单位：mm)；Rs_i为时段内地表径流量(单位：mm)；ETa_i为时段内蒸散发量(单位：mm)；WUBi为时段内表层土壤向深层的下渗量(单位：mm)。
对于深层土壤，在接纳上层土壤的下渗水分的同时，产生深层壤中流，同时植被根系吸水用于蒸腾作用，深层土壤的水量平衡方程如下：
WUBi+AWd_i＝AWd_i+1+Rsd_i+ETad_i+WUD_i
式中：AWd_i为时段初深层土壤的含水量(单位：mm)；AWd_i+1为时段末深层土壤的含水量(单位：mm)；Rsd_i为深层土壤出流量(单位：mm)；ETad_i为时段内植被蒸腾所耗的深层土壤水分(单位：mm)；WUD_i为土壤向地下的入渗量(单位：mm)。
由此，以上水量平衡方程可表达为土壤含水量的函数，对水量平衡方程中的土壤含水量求导，通过牛顿跌代公式计算时段内的土壤含水量：
x1=x0-f(x0)f′(x0)]]>
求算出土壤含水量，即可计算时段内各水文分量，包括蒸散发、地表径流、表层壤中流、深层壤中流和地下径流。
步骤S2，建立非点源污染模型，对农业区非点源污染进行诊断。
农业区污染包括点源和非点源污染，其中工厂和企业形成的点源污染可以通过统计数据获得，而由农业种植、养殖、居民点等形成的非点源污染是水环境污染诊断的核心模型：
非点源污染总量NP可以表示为：
NP＝P_ol-P_oi，
式中，P_ol为农业区(或子区)污染总量，P_oi为点源污染负荷，NP为非点源污染负荷。所以，子区域水环境诊断值S_i可表达为：
S_i＝NP_i/R_i，式中，NP_i为非点源污染量，R_i为产生的径流量。
根据S_i值诊断是否非点源污染物超排放标准。
对于非点源污染的计算，以下通过实施例具体描述。
首先进行污染源强的计算。
污染源强即为各污染类型的排放量，是溶解态污染负荷计算的一个参数。本文根据郑东海等的研究成果进行农田生产、农村居民点、城市径流和畜禽养殖4种污染类型的污染源强计算。
S71，计算农业生产的污染源强：
Qa=ΣxnΣymAx&CenterDot;ωxy&CenterDot;M]]>
其中，x为农田类型，n为农田类型数，y为农田污染物指标类型，主要包括总氮、总磷和氨氮3种类型，m为农田染物指标类型的类型数，A为农田面积，ω_xy为各种作为不同污染物指标的源强系数，M为修正系数，包括坡度、土壤、化肥和降水量修正；
S72，计算畜禽养殖污染源强：
Qb=ΣhnnΣkmmCh&CenterDot;Dh&CenterDot;kh&CenterDot;ωhk&CenterDot;ηhk]]>
其中，h为畜禽养殖类型，主要包括大牲畜、猪、羊和家禽4种类 型；nn为畜禽养殖类型数，k为畜禽养殖污染物指标类型，mm为畜禽养殖污染物指标类型数，C为畜禽数量，D为不同畜禽类型饲养周期，k_h为不同畜禽类型粪便排泄指数，ω_hk为不同畜禽类型不同污染物指标的源强系数；η_hk为不同畜禽类型不同污染物指标的流失率。
S73，计算农村居民点的污染源强：
Qc=ΣfmnPc&CenterDot;ωf,]]>
其中，f为污染物指标类型，主要包括总氮、总磷、COD和氨氮4种类型，mn为类型数，P_c为乡村人口数量，ω_f为不同污染物指标的源强系数。
S74，计算城镇居民点的污染源强：
Qd=ΣgnmPd&CenterDot;ωg,]]>
其中，g为污染物指标类型，nm为类型数，P_d为城镇人口数量，ω_g为不同污染物指标的源强系数。
S75，计算所述农业区的污染源强：Q＝Qa+Qb+Qc+Qd。
非点源污染的计算，具体包括以下步骤：
S21，计算溶解态非点源污染的值。具体公式为：
C=Σi=14&epsiv;&epsiv;0&CenterDot;(1-ekR0t)&CenterDot;Qi&CenterDot;Ni&CenterDot;Si,]]>
其中，C为污染负荷，i为溶解态污染类型，n为类型数；Q_i为污染源强；k为地面冲刷系数；R₀为标准雨强,此处可以取值12.7；ε为径流系数，ε₀为标准径流系数，反映不透水硬化地面，此处可以取值0.87；N_i为自然因子修正系数，如坡度、植被覆盖等因子；S_i为社会因子修正系数，反映社会发展程度对非点源污染源强的削弱程度。
其中自然修正因子可以表示为：
N＝slope_CO·vegcover_CO
slopeCO=slope-slopeminslopemax-slopemin]]>
veg cov erCO=veg cov ermax-veg cov erveg cov ermax-veg cov ermin]]>
式中，slope_CO和vegcover_CO分别为坡度和植被覆盖度的修正系数；slope_max、slope_min和slope分别为最大坡度、最小坡度和计算单元坡度；vegcover_max、vegcover_min和vegcover分别为最大植被覆盖度、最小植被覆盖度和计算单元植被覆盖度。
社会修正因子可以表示为：
S＝(1-W)·(1-U)
式中，W和U分别表示垃圾处理率和非点源污染进入管网比率。
S22，计算吸附态非点源污染的值：
基于通用土壤流失方程，建立了吸附态氮和磷流失量的估算模型如下：
C_a＝X·Q_a·τ
式中,X为土壤侵蚀量；Q_a为流失土壤中的氮磷污染物背景含量，通过田间调查获取的不同土壤类型氮磷含量计算；τ为吸附态污染物富集系数；
其中τ为吸附态污染物富集系数的计算公式为：
lnτ＝2-0.2lnX。
在上述的计算中，溶解态二元结构非点源污染模型及吸附态非点源污染模型计算过程中用到数据主要包括社会经济数据(乡村人口数量、城镇人口数量和畜禽数量)、社会参数(垃圾处理率和非点源污染进入管网率)、气象数据(降雨时长)、环境参数(土地利用/覆盖分类图、坡度、植被覆盖度、土壤侵蚀量)等。
(1)社会经济数据和参数
农村居民点、城市地表径流和畜禽养殖是主要的非点源污染源。乡村人口数量、城镇人口数量和畜禽养殖规模直接控制着非点源污 染源强和规模，获取准确的污染源强数据是进行非点源污染模拟的前提。乡村人口数量、城镇人口数量和畜禽养殖数量主要依据全国和省市统计年鉴。
社会参数主要是垃圾处理率和非点源污染进入管网率。由于不同污染源所处的地区基础设施建设情况的不同，直接导致了非点源污染物流失量的百分比不同。我国城市地区基础设施建设相对完善，部分垃圾要回收处理，大部分污水排往污水管道；而我国农村的基础设施建设相对薄弱，缺乏排水管网和垃圾回收处理设施。
(2)气象数据和坡度
非点源污染物是伴随着地表径流迁移转化的，而降雨事件是地表径流的主要驱动因子，降雨数据是非点源污染负荷模拟的主要参数之一。目前降雨数据的获取方式多种多样，包括地面气象站点实时监测以及通过气象卫星反演等方式。
地形坡度直接影响着土壤侵蚀和地表径流，是非点源污染模拟的重要自然因子。目前主要通过DEM数据或地形图等提取地形坡度参数。
(3)植被覆盖度
地表植被有保水保土保肥的功效，可以拦截雨滴、改良土壤抗蚀性能，大大减少降雨或融雪径流的产流和产沙量，拦截污染物质随径流继续迁移，从而有效控制非点源污染的产生。地表植被的覆盖常常用植被覆盖度来衡量，植被覆盖度是溶解态二元结构非点源污染模型和土壤侵蚀量计算的重要参数。植被覆盖度的遥感分析方法主要有混合像元分解法和经验模型法。
(4)土地利用/覆盖分类
土地利用/覆盖数据控制着非点源污染源和污染物种类的空间分布情况，也影响污染物质、土壤和径流的输入、输出和迁移的过程。同时，土地利用/覆盖状况是人工干预最易实现的环节，也是水 环境保护的重要手段，因此土地利用/覆盖是非点源污染研究的重要参数。目前基于遥感技术的土地利用/覆盖分类技术比较成熟，主要有基于传统统计分析的分类方法(监督分类和非监督分类)、神经网络分类方法、模糊分类方法、决策树分类方法和专家系统分类方法等。
(5)土壤侵蚀量
土壤侵蚀是破坏农田设施，降低土壤肥力，淤积水库、渠道，抬高河床等的主要元凶之一，同时，土壤营养物质吸附在土壤颗粒伴随土壤侵蚀流失，富集在河流、水库和湖泊等水体中，是非点源污染的主要形式之一。土壤侵蚀量是计算吸附态非点源污染的主要参数之一。目前土壤侵蚀研究方法较多，基于遥感技术的土壤侵蚀研究成为目前的热点，关于全国、重点区域和流域的土壤侵蚀遥感调查成果和数据库建设的报道屡见不鲜。
以上参数的获取方案如下：
1)社会经济数据(乡村人口数量、城镇人口数量和畜禽养殖数量)主要通过《中国统计年鉴·2011》《河北省统计年鉴·2011》《山西省统计年鉴·2011》和《河南省统计年鉴·2011》获取；
2)社会参数(垃圾处理率和非点源污染进入管网率)主要通过相关文献研究和统计年鉴获取；
3)气象数据主要通过研究区各县市气象站地面监测获取；
4)坡度数据主要是通过DEM数据利用ARCGIS软件提取；
5)土地利用/覆盖分类、植被覆盖度、土壤侵蚀量主要利用遥感数据、DEM数据、气象数据、土壤数据等进行遥感反演和模拟。
S23，计算非点源污染量的值：
NP＝C+C_a。
步骤S3，根据所述农业区用水的诊断和所述非点源污染的诊断，对农业区水资源环境进行诊断。
本发明的另一个实施例中，提供一种农业区水资源环境诊断系统。
图2示出了本发明的一种农业区水资源环境诊断系统的结构框图。
参照图2，本发明的实施例的农业区水资源环境诊断系统包括：
农业区用水诊断模块10，用于建立蒸散发模型和产流模型，根据所述蒸散发模型和所述产流模型对农业区用水进行诊断；
农业区非点源污染诊断模块20，用于建立非点源污染模型，对农业区非点源污染进行诊断；
农业区水资源环境诊断模块30，根据所述农业区用水的诊断和所述非点源污染的诊断，对农业区水资源环境进行诊断。
通过上述实施例可知，本发明的一种农业区水资源环境的诊断方法及其系统，基于蒸散发模型、径流模型和非点源污染模型，以年为时间尺度，并综合了径流数据，有效地解决了农业区不同种植结构总用水和灌溉用水估算问题和非点源污染输出强度估算问题。同时解决了传统非遥感技术估算方法中因未考虑下垫面差异的问题，并最终能够有效地诊断水资源环境问题。
虽然结合附图描述了本发明的实施方式，但是本领域技术人员可以在不脱离本发明的精神和范围的情况下做出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。