技术领域
[0001] 本
发明涉及
灌溉领域，特别涉及一种灌水频率与强度的确定方法。
背景技术
[0002] 农田灌水频率指灌水间隔的时长，农田灌水强度指每次灌水的水量，农田灌水频率和灌水强度是拟定农田灌溉计划需要确定的重要内容。
[0003] 少数有条件的地区，通常采用农田灌溉试验方法来确定灌水频率和灌水强度，其他地区只能依靠以往的灌水经验确定。然而，试验方法确定灌水频率和灌水强度是一项费时费
力的工作，不同地区的灌水频率和灌水强度也不能通用，因此，目前大部分地区仅凭经验来拟定农田灌溉计划，造成水资源利用效率低下。
发明内容
[0004] 为了提高水资源的利用效率，本发明提供了一种灌水频率与强度的确定方法。所述技术方案如下：
[0005] 本发明提供了一种灌水频率与强度的确定方法，所述方法包括：
[0006] 接收用户输入的多个灌水强度集合、下边界条件方程的方程名称、T个单位时段内的每个单位时段的理论蒸腾强度、所述每个单位时段的下边界水分通量和待测试
土壤中的Z个
位置点中每个位置点在第0单位时段的土壤含水率，所述Z个位置点位于同一条直线上且所述直线垂直于所述待测试土壤的表面，所述Z个位置点中相邻两位置点之间的间隔为预设距离，所述多个灌水强度集合中的任一个灌水强度集合包括多个灌水强度以及所述多个灌水强度中的每个灌水强度对应的灌水持续时间和灌水起始时间；
[0007] 建立所述方程名称对应的下边界条件方程；
[0008] 根据目标灌水强度集合中包括每个灌水强度对应的灌水起始时间，确定所 述目标灌水强度集合中是否定义在第i个单位时段内灌水，所述目标灌水强度集合是所述多个灌水强度集合中的任一个灌水强度集合，i＝1、2、3……T，T为大于0的整数；
[0009] 如果没有定义在第i个单位时段内灌水，则根据所述每个位置点在第i-1单 位时段的土壤含水率计算第i个单位时段内待测试土壤的表层平均体积含水率 根据所述表层平均体积含水率 构建上边界方程，且所述上边界方程为如下公式(1)所示的第i单位时段的目标上边界方程：
[0010]
[0011] 在公式(1)中，ES(i)为第i单位时段内的土壤
蒸发强度；θ1和θ2分别为预 设的表层土壤平均体积含水率，ES0为预设水面蒸发强度，a、b为预设值；
[0012] 如果定义在第i个单位时段内灌水，根据所述每个位置点在第i-1单位时段 的土壤含水率建立第i单位时段的目标上边界方程；
[0013] 根据所述下界条件方程、所述目标上边界方程和预设的土壤含水率分布变 化模型计算第i个单位时段内所述每个位置点的土壤含水率；
[0014] 根据第i个单位时段内所述每个位置点的土壤含水率，按如下公式(2)计 算第i个单位时段内作物的实际蒸腾量
[0015]
[0016] 在公式(2)中， 为第z个位置点在第i个单位时段内的土壤含水率，L(i) 为作物根系在第i个单位时段内的生长长度，EC(i)为第i单位时段的理论蒸腾强 度，Δz为相邻两个位置点之间的距离，Δt为单位时段的时长，A、θmax和θmax为 别为预设值；
[0017] 根据T个单位时段中的每个单位时段内作物的实际蒸腾量 和所述目标 灌水强度集合，按如下公式(3)计算所述目标灌水强度集合对应的水分利用率 fm；
[0018]
[0019] 在公式(3)中，I(i)是所述目标灌水强度集合中定义的第i个单位时段内 的灌水强度，t(i)为第i单位时段内的灌水持续时间；
[0020] 从所述多个灌水强度集合中选择最大水分利用率对应的灌水强度集合；
[0021] 输出所述选择的灌水强度集合包括的每个灌水强度以及所述每个灌水强度 对应的灌水持续时间和灌水起始时间。
[0022] 可选的，所述建立所述方程名称对应的下边界条件方程，包括：
[0023] 在所述方程名称为定
水头边界时，建立如下公式(4)所示的第一下边界条 件方程为：
[0024]
[0025] 在公式(4)式中，hz为下边界定水头，为预设值， 为第Z个位置点在第 i单位时段的土壤
负压水头；
[0026] 在所述方程名称为通量边界，建立如下公式(5)所示的第二下边界条件方 程为：
[0027]
[0028] 在公式(5)中，h为土壤负压水头，α为所述待测试土壤与水平面之间的 夹
角；K为土壤导水率；σ(i)为第i单位时段内的下边界水分通量；其中， Ks为
土壤饱和导水率，为预设值，h为土壤负压水头；θs、θr分别为饱和含水率 和残余含水率，均为预设值；a、n均为预设参数；
[0029] 在所述方程名称为自由排水边界，建立如下公式(6)所示的第三下边界条件方程为：
[0030]
[0031] 可选的，所述根据所述每个位置点在第i-1单位时段的土壤含水率计算第i个单位时段内待测试土壤的表层平均体积含水率 根据所述表层平均体积含水率 构建上边界方程，且所述上边界方程为如下公式(1)所示的第i单位时段的目标上边界方程，包括：
[0032] 根据第1个位置点在第i-1单位时段的土壤含水率、第2个位置点在第i-1单位时段的土壤含水率和第3个位置点在第i-1单位时段的土壤含水率，计算出第i个单位时段内待测试土壤的表层平均体积含水率 并建立如公式(1)所示的第一上边界条件方程，其中第1个位置点位于所述待测试土壤的表面；
[0033] 根据所述下边界条件方程、如公式(1)所示的第一上边界条件方程和如下公式(7)所示的预设的土壤含水率分布变化模型，计算出所述每个位置点在第i单位时段的中间土壤含水率
[0034]
[0035] 根据第1个位置点在第i单位时段的中间土壤含水率、第2个位置点在第i单位时段的中间土壤含水率和第3个位置点在第i单位时段的中间土壤含水率，重新计算出第i个单位时段内待测试土壤的表层平均体积含水率 根据重新计算出的表层平均体积含水率重新构建如公式(1)所示的目标上边界方程。
[0036] 可选的，所述如果定义在第i个单位时段内灌水，根据所述每个位置点在第i-1单位时段的土壤含水率建立第i单位时段的目标上边界方程，包括：
[0037] 建立如下公式(8)所述的第二上边条件方程，
[0038]
[0039] 根据所述下边界条件方程、如公式(8)所示的第二上边界条件方程和如公式(7)所示的预设的土壤含水率分布变化模型，计算出所述每个位置点在第i单位时段的中间土壤含水率
[0040] 如果第1个位置点在第i单位时间段的的中间土壤含水率 小于预设值θs，则确定地表土壤含水率是未达到饱和，则将所述公式(8)所述的第二上边条件方程确定为目标上边条件方程。
[0041] 可选的，所述方法还包括：
[0042] 如果地表土壤含水率达到饱和，则构建的目标上边界条件方程为如下公式 (9)所示的第三上边界方程：
[0043]
[0044] 式中， 为第1个位置点在第i+1个单位时段内的水头；k1、k2为地表积水开始和结束时刻；zxjl(i+1)、zxjl(i)分别为第i+1和第i单位时段内地表的积水深度；ZXNL为地表植被滞蓄能力，为预设数值；RS(i+1)为第i+1单位时段内的入渗强度，dt(i+1)为第i+1单位时段的时间步长；当地表积水深
度zxjl(i+1)小于地表植被滞蓄能力ZXNL时，受植被滞蓄作用新产生的超渗产流滞蓄在地表，形成静水压力水头，数值上与积水深度相等；当地表积水深度zxjl(i+1)大于地表植被滞蓄能力ZXNL时，超过地表滞蓄能力的积水形成径流流走，此时地表积水深度达到最大，数值上与地表植被滞蓄能力相等。
[0045] 可选的，所述根据所述下界条件方程、所述目标上边界方程和预设的土壤含水率分布变化模型计算第i个单位时段内所述每个位置点的土壤含水率，包括：
[0046] 第一步：将第i-1个单位时段内所述每个位置点的土壤含水率作为第i个单位时段内所述每个位置点的第一土壤含水率
[0047] 第二步：根据所述下界条件方程、所述目标上边界方程和预设的土壤含水率分布变化模型计算第i个单位时段内所述每个位置点的土壤含水率 将所述第i个单位时段内所述每个位置点的土壤含水率 作为第二土壤含水率
[0048] 第三步：确定第i个单位时段内所述每个位置点的第一土壤含水率 和第二土壤含水率 是满足如下公式(10)所示的条件，ε为允许的相对误差，为预设数值；
[0049]
[0050] 如果满足公式(10)所示的条件，则确定每个位置点在第i个单位时段内的土壤含水率 等于第二含水率 即
[0051] 如果不满足公式(10)所示的条件，则将每个位置点在第i个单位时段内的第二含水率 分别作为在第i个单位时段内的第一含水率 即 然后返
回执行第二步。
[0052] 本发明提供的技术方案的有益效果是：
[0053] 在本发明
实施例中，通过接收用户输入的多个灌水强度集合，然后计算出每个灌水强度集合对应的水分利用率，水分利用率越高的灌水强度集合，造成的水资源利用率越高，水分利用率越低的灌水强度集合，造成的水资源利用率越低。在本实施例中，输出水分利用率最高的灌水强度集合包括的每个灌水强度以及每个灌水强度对应的灌水持续时间和灌水起始时间，供用户对作物进行灌水，从而提高了水资源的利用效率。
附图说明
[0054] 图1是本发明实施列提供了一种灌水频率与强度的确定方法
流程图；
[0055] 图2是本发明实施例提供的土壤剖面含水率分布图；
[0056] 图3是本发明实施例提供的土壤剖面负压水头分布图。
具体实施方式
[0057] 为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
[0058] 实施例
[0059] 参见图1，本发明实施列提供了一种灌水频率与强度的确定方法，所述方法包括：
[0060] 步骤101：接收用户输入的多个灌水强度集合、下边界条件方程的方程名称、T个单位时段内的每个单位时段的理论蒸腾强度、每个单位时段的下边界水分通量和待测试土壤中的Z个位置点中每个位置点在第0单位时段的土壤含水率，Z个位置点位于同一条直线上且直线垂直于所述待测试土壤的表面，Z个位置点中相邻两位置点之间的间隔为预设距离，多个灌水强度集合中的任一个灌水强度集合包括多个灌水强度以及多个灌水强度中的每个灌水强度对应的灌水持续时间和灌水起始时间。
[0061] 本实施例的执行主体可以为计算机，用户可以通过计算机的
键盘等输入输出设备向计算机输入多个灌水强度集合、初始条件。
[0062] 用户除了向计算机中输入上述几个参数外，还可以输入蒸发参数、空间步长等其他参数。例如，参见表1，用户可以输入的参数可以为。
[0063] 表1
[0064]
[0065]
[0066] 步骤102：建立所述方程名称对应的下边界条件方程。
[0067] 下边界条件方程有第一下边界条件方程、第二下边界条件方程和第三下边界条件方程三种，用户根据实际情况选择一种下边界条件方程并向计算机输入该下边界条件方程的方程名称。
[0068] 本步骤可以通过如下三个步骤来实现，分别为：
[0069] 1021：在该方程名称为定水头边界时，建立如下公式(4)所示的第一下边界条件方程为：
[0070]
[0071] 在公式(4)式中，hz为下边界定水头，为预设值，为第Z个位置点在第i单位时段的土壤负压水头；
[0072] 1022：在该方程名称为通量边界，建立如下公式(5)所示的第二下边界条件方程为：
[0073]
[0074] 在公式(5)中，h为土壤负压水头，α为待测试土壤与水平方向的夹角；K为土壤导水率；σ(i)为第i单位时段内的下边界水分通量，第i个单位时段内的 下边界水分
能量σ(i)是用户输入事先输入的；其中，Ks为土壤饱和导水率，为预设值，h为土壤负压水头；θs、θr分别为饱和含水率 和残余含水率，均为预设值；a、n均为预设参数；
[0075] 1023：在该方程名称为自由排水边界，建立如下公式(6)所示的第三下边 界条件方程为：
[0076]
[0077] 步骤103：根据目标灌水强度集合中包括每个灌水强度对应的灌水起始时 间，确定目标灌水强度集合中是否定义在第i个单位时段内灌水，目标灌水强度 集合是多个灌水强度集合中的任一个灌水强度集合，i＝1、2、3……T，T为大于 0的整数。
[0078] 步骤104：如果没有定义在第i个单位时段内灌水，则根据每个位置点在第 i-1单位时段的土壤含水率计算第i个单位时段内待测试土壤的表层平均体积含 水率 根据该表层平均体积含水率 构建上边界方程，且该上边界方程为如下 公式(1)所示的第i单位时段的目标上边界方程。
[0079]
[0080] 在公式(1)中，ES(i)为第i单位时段内的土壤蒸发强度；θ1和θ2分别为预 设的表层土壤平均体积含水率，ES0为预设水面蒸发强度，a、b为预设值。
[0081] 本步骤可以为：①，首先根据第1个位置点在第i-1单位时段的土壤含水率、 第2个位置点在第i-1单位时段的土壤含水率和第3个位置点在第i-1单位时段 的土壤含水率，计算出第i个单位时段内待测试土壤的表层平均体积含水率 并根据该表层平均体积含水率先初步建立如公式(1)所示的第一上边界条件方程，其中第1个位置点位于待测试土壤的表面。
[0082] 需要说明的是：此时由于该表层平均体积含水率 是根据第i-1个单位时段 内的土壤含水率计算得到的，因此计算出的该表层平均体积含水率 可能不是第 i单位时段内真实值，导致初步建立的如公式(1)所示的第一上边界条件方程 可能有误。为此，本步骤还需要通过如下方式进行校正，校正过程如下：
[0083] ②，然后根据该下边界条件方程、如公式(1)所示的第一上边界条件方程 和如下公式(7)所示的预设的土壤含水率分布变化模型，计算出每个位置点在 第i单位时段的中间土壤含水率
[0084]
[0085] A为预设值。
[0086] 例如，假设在步骤102中建立的下边界条件方程为如公式(5)所示的第二下边界条件方程，则根据公式(1)、(5)和(7)计算出每个位置点在第i单位时段的中间土壤含水率[0087] ③，最后根据第1个位置点在第i单位时段的中间土壤含水率、第2个位置点在第i单位时段的中间土壤含水率和第3个位置点在第i单位时段的中间土壤含水率，重新计算出第i个单位时段内待测试土壤的表层平均体积含水率 根据重新计算出的表层平均体积含水率 重新构建如公式(1)所示的第一上边界条件方程，并作为目标上边界方程。
[0088] 步骤105：如果定义在第i个单位时段内灌水，根据每个位置点在第i-1单位时段的土壤含水率建立第i单位时段的目标上边界方程。
[0089] 在第i个单位时段内灌水包括土壤表层的土壤含水率饱和和未饱和两种情况，在每种情况下的建立的目标上边界条件的方程不同，详细建立过程如下：
[0090] 1051：先假设土壤表层的土壤含水率未饱和，即第1个位置点的未达到饱和，在未饱和情况下建立如下公式(8)所述的第二上边条件方程，
[0091]
[0092] 1052：根据该下边界条件方程、如公式(8)所示的第二上边界条件方程和如公式(7)所示的预设的土壤含水率分布变化模型，计算出每个位置点在第i单位时段的中间土壤含水率
[0093] 例如，假设在步骤102中建立的下边界条件方程为如公式(5)所示的第二下边界条件方程，则根据公式(5)、(8)和(7)计算出每个位置点在第i单位时段的中间土壤含水率[0094] 1053：如果第1个位置点在第i单位时间段的的中间土壤含水率 小于预设值θs，则确定地表土壤含水率是未达到饱和，则将公式(8)所述的第二上边条件方程确定为目标上边条件方程。
[0095] 1054：如果地表土壤含水率达到饱和，则构建的目标上边界条件方程为如下公式(9)所示的第三上边界方程：
[0096]
[0097] 式中， 为第1个位置点在第i+1个单位时段内的水头；k1、k2为地表积水开始和结束时刻；zxjl(i+1)、zxjl(i)分别为第i+1和第i单位时段内地表的积水深度；ZXNL为地表植被滞蓄能力，为预设数值；RS(i+1)为第i+1单位时段内的入渗强度，dt(i+1)为第i+1单位时段的时间步长；当地表积水深
度zxjl(i+1)小于地表植被滞蓄能力ZXNL时，受植被滞蓄作用新产生的超渗产流滞蓄在地表，形成静水压力水头，数值上与积水深度相等；当地表积水深度zxjl(i+1)大于地表植被滞蓄能力ZXNL时，超过地表滞蓄能力的积水形成径流流走，此时地表积水深度达到最大，数值上与地表植被滞蓄能力相等。
[0098] 步骤106：根据该下界条件方程、该目标上边界方程和公式(7)所示的预设的土壤含水率分布变化模型计算第i个单位时段内所述每个位置点的土壤含水率。
[0099] 本步骤可以为：
[0100] 第一步：将第i-1个单位时段内每个位置点的土壤含水率作为第i个单位时段内所述每个位置点的第一土壤含水率
[0101] 第二步：根据该下界条件方程、该目标上边界方程和如公式(7)所示的预设的土壤含水率分布变化模型计算第i个单位时段内所述每个位置点的土壤含水率 将第i个单位时段内所述每个位置点的土壤含水率 作为第二土壤含水率
[0102] 例如，假设在步骤102中建立的下边界条件方程为如公式(5)所示的第二下边界条件方程，在步骤105中建立的目标上边界条件方程是如公式(8)所示的第二上边界条件方程，则根据公式(5)、(8)和(7)计算出每个位置点在第i单位时段的土壤含水率 作为第二土壤含水率
[0103] 第三步：确定第i个单位时段内每个位置点的第一土壤含水率 和第二土壤含水率 是满足如下公式(10)所示的条件，ε为允许的相对误差，为预设数值；
[0104]
[0105] 如果满足公式(10)所示的条件，则确定每个位置点在第i个单位时段内的土壤含水率 等于第二含水率 即 结束返回；
[0106] 如果不满足公式(10)所示的条件，则将每个位置点在第i个单位时段内的第二含水率 分别作为在第i个单位时段内的第一含水率 即 然后返回执行第二步。
[0107] 步骤107：根据第i个单位时段内每个位置点的土壤含水率，按如下公式(2)计算第i个单位时段内作物的实际蒸腾量
[0108]
[0109] 在公式(2)中， 为第z个位置点在第i个单位时段内的土壤含水率，L(i) 为作物根系在第i单位时段内的生长长度，EC(i)为第i单位时段的理论蒸腾强度， Δz为相邻两个位置点之间的距离，Δt为单位时段的时长，θmax和θmax为别为预设 值。其中，L(i)＝-4.40132
+0.7744i+0.00036i。
[0110] 重复执行上述步骤103至107的步骤计算出第1个单位时段的实际蒸腾量 第2个单位时段的实际蒸腾量 直至计算出第T个单位时段的实际 蒸腾量 时为
止。
[0111] 步骤108：根据T个单位时段中的每个单位时段内作物的实际蒸腾量 和 目标灌水强度集合，按如下公式(3)计算目标灌水强度集合对应的水分利用率 fm；
[0112]
[0113] 在公式(3)中，I(i)是目标灌水强度集合中定义的第i个单位时段内的灌 水强度，t(i)为第i单位时段内的灌水持续时间。
[0114] 对于其他每个灌水强度集合，同上述目标灌水强度集合一样按上述步骤103 至108的操作计算出其他每个灌水强度集合对应的水分利用率fm。
[0115] 步骤109：从该多个灌水强度集合中选择最大水分利用率对应的灌水强度集 合。
[0116] 步骤110：输出选择的灌水强度集合包括的每个灌水强度以及每个灌水强度 对应的灌水持续时间和灌水起始时间。
[0117] 在本发明实施例中，通过接收用户输入的多个灌水强度集合，然后计算出 每个灌水强度集合对应的水分利用率，水分利用率越高的灌水强度集合，造成 的水资源利用率越高，水分利用率越低的灌水强度集合，造成的水资源利用率 越低。在本实施例中，输出水分利用率最高的灌水强度集合包括的每个灌水强 度以及每个灌水强度对应的灌水持续时间和灌水起始时间，供用户对作物进行 灌水，从而提高了水资源的利用效率。
[0118] 如下本发明提供的应用实例。
[0119] 在该应用实例中，以华北某小麦试验田为例，允许最大含水率、最小含水 率分别取0.296、0.08，土壤水动力参数见下表2。
[0120] 表2
[0121]
[0122]
[0123] 考虑根系吸水速率的
时空分布，根系吸水率与吸水深度呈线性变化，计算 公式如下：
[0124]
[0125] 式中，EC(t)为供水充分条件下植株蒸腾强度；A为经济系数，取0.597；t 为自冬小麦
播种日起至计算起始时间为止的天数。当t取较小时zr(t)<0，与实际 不符，取zr(t)＝0。
[0126] 棵间蒸发强度考虑表土5cm内土壤平均含水率对蒸发的影响，由冬小麦田 间试验得出以下关系：
[0127]
[0128] 式中，为
表土层5cm内平均体积含水率；ES0为水面蒸发强度。
[0129] 充分供水条件冬小麦理论蒸腾强度的基本资料，见下表3：
[0130] 表3
[0131]生育期 播种-越冬 越冬-返青 返青-拔节 拔节-抽穗 抽穗-成熟
天数/d 78 114 153 185 230
ES0/mm/d 0.53 0.28 0.53 0.27 0.23
EC/mm/d 0.37 0.57 1.23 3.52 3.78
[0132] 空间步长取4cm，模拟土层深度位200cm，下边界条件为自由排水边界， 模拟时长为整个小麦生育230d，时间步长取0.01d，地表最大积水深度取1cm， 地表坡度为0，最大允许误差取0.0001。
[0133] 接收输入的灌水频率与强度，即得到如下表4所示的灌水强度集合：
[0134] 表4
[0135]
[0136] 下表5为计算得到的各个组合的水分利用率，由此可知灌水频率为30d，灌 水强度为4.4cm的组合最优。
[0137] 表5
[0138]灌水强度集合 水分利用率
集合1 0.6514
集合2 0.6596
集合3 0.6935
集合4 0.8216
[0139] 图2和图3为最优灌水频率和强度下的第114天、153天、230天土壤剖面 含水率分布和土壤负压水头分布。
[0140] 本领域普通技术人员可以理解实现上述实施例的全部或部分步骤可以通过
硬件来完成，也可以通过程序来指令相关的硬件完成，所述的程序可以存储于 一种计算机可读存储介质中，上述提到的存储介质可以是只读
存储器，磁盘或 光盘等。
[0141] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的 精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的 保护范围之内。
