CN113255093B - 无资料地区径流估算方法、装置、计算机设备和存储介质 - Google Patents

无资料地区径流估算方法、装置、计算机设备和存储介质 Download PDF

Info

Publication number
CN113255093B
CN113255093B CN202110065530.4A CN202110065530A CN113255093B CN 113255093 B CN113255093 B CN 113255093B CN 202110065530 A CN202110065530 A CN 202110065530A CN 113255093 B CN113255093 B CN 113255093B
Authority
CN
China
Prior art keywords
parameter
value
function
water
objective function
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110065530.4A
Other languages
English (en)
Other versions
CN113255093A (zh
Inventor
李雪莹
龙笛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Original Assignee
Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tsinghua University filed Critical Tsinghua University
Priority to CN202110065530.4A priority Critical patent/CN113255093B/zh
Publication of CN113255093A publication Critical patent/CN113255093A/zh
Application granted granted Critical
Publication of CN113255093B publication Critical patent/CN113255093B/zh
Priority to US17/697,915 priority patent/US20230030888A1/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N33/00—Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
    • G01N33/24—Earth materials
    • G01N33/246—Earth materials for water content
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/10—Devices for predicting weather conditions

Landscapes

  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Chemical & Material Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Geology (AREA)
  • Atmospheric Sciences (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Ecology (AREA)
  • Environmental Sciences (AREA)
  • Pathology (AREA)
  • Remote Sensing (AREA)
  • General Health & Medical Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Food Science & Technology (AREA)
  • Medicinal Chemistry (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)

Abstract

本申请涉及一种无资料地区径流估算方法、装置、计算机设备和存储介质。所述方法包括：通过获取驱动数据，所述驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；然后，构建流域月尺度水量平衡变化的目标函数，并确定目标函数的待求解参数；对待求解参数进行初始化，将待求解参数的初始值与驱动数据输入目标函数，得到目标函数的初始函数值；并基于所述待求解参数的初始值和所述初始函数值，对所述待求解参数进行迭代计算；最后，根据所述待求解参数的优化值估算预设时间段内流域月尺度径流值。采用本方法能够在不依赖任何实测径流资料的条件下，通过输入较少的驱动数据即可对无资料地区的径流值进行估算，估算准确率高。

Description

无资料地区径流估算方法、装置、计算机设备和存储介质
技术领域
本申请涉及水文水资源管理技术领域，特别是涉及一种无资料地区径流估算方法、装置、计算机设备和存储介质。
背景技术
河流径流量是水文循环中的重要变量，也是用于区域水资源评估、生态环境保护和气候变化分析的基础数据。相关的径流监测数据主要来源于水文监测站点，虽然其精度较高，但是全球站点分布十分不均匀且多数站点观测时段较短，受地形及观测成本的制约，许多无资料区缺少可靠的径流数据。
目前的研究通常采用搭建水文模型的方法模拟流域径流量，然而相关方法严重依赖实测径流资料对模型参数进行率定，在无资料地区应用受限。
目前使用相邻流域的实测流量数据率定该流域的模型参数并移植至无资料地区，但是由于不同流域的独特性，模拟准确率不高。
发明内容
基于此，有必要针对上述技术问题，提供一种能够提高无资料地区径流估算准确率的无资料地区径流估算方法、装置、计算机设备和存储介质。
第一方面，提供了一种无资料地区径流估算方法，该方法包括：
获取驱动数据，该驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；
构建流域月尺度水量平衡变化的目标函数，确定该目标函数的待求解参数；
对该待求解参数进行初始化，得到待求解参数的初始值；
将该待求解参数的初始值和该驱动数据输入该目标函数，得到目标函数的初始函数值；
基于该待求解参数的初始值和该初始函数值，对该待求解参数进行迭代计算，得到待求解参数的优化值，将该待求解参数的优化值输入该目标函数得到目标函数值；
在该目标函数值满足预设的迭代终止条件的情况下，根据该待求解参数的优化值估算预设时间段内流域月尺度径流值。
在其中一个实施例中，对该待求解参数进行迭代计算，得到待求解参数的优化值，包括：
计算该目标函数对该待求解参数的梯度函数；
计算该梯度函数在上一轮待求解参数的优化值下的梯度值；
计算本轮该待求解参数在该梯度值下的参数值，得到该待求解参数的优化值。
在其中一个实施例中，计算本轮该待求解参数在该梯度值下的参数值，得到该待求解参数的优化值，包括：
根据该上一轮待求解参数的优化值、该梯度值以及预设的优化精度，将所述上一轮待求解参数的优化值沿着所述梯度值的负方向进行迭代计算，得到该待求解参数的优化值。
在其中一个实施例中，该待求解参数包括多个待求解参数，共同组成参数空间；
计算该目标函数对该待求解参数的梯度函数，包括：
计算该目标函数对该参数空间中每一个待求解参数的梯度函数；
计算该梯度函数在上一轮待求解参数的优化值下的梯度值，包括：
计算该每一个待求解参数的梯度函数在对应的上一轮待求解参数的优化值下的梯度值；
根据该上一轮待求解参数的优化值、该梯度值以及优化精度，将该上一轮待求解参数的优化值沿着该梯度值的负方向进行迭代计算，得到该待求解参数的优化值，包括：
根据该上一轮每一个待求解参数的优化值、对应的该梯度值以及优化精度，将该上一轮每一个待求解参数的优化值沿着该梯度值的负方向进行迭代计算，得到对应的优化后的参数值，该每一个优化后的参数值共同组成优化的参数空间；
所述将所述待求解参数的优化值输入所述目标函数得到目标函数值，包括：
将所述优化的参数空间输入所述目标函数得到目标函数值。
在其中一个实施例中，构建流域月尺度水量平衡变化的目标函数，包括：
根据水量平衡模型中水量的通量变化关系进行转换，得到第一土壤水储量变化参数；
根据该水量平衡模型中水量的状态量变化关系，得到第二土壤水储量变化参数；
根据该第一土壤水储量变化参数、该第二土壤水储量变化参数构建该流域月尺度水量平衡变化的目标函数。
在其中一个实施例中，该水量平衡模型的构建过程包括：
根据预设时间段内水量输入通量函数、水量输出通量函数以及水量的状态量之间的关系构建该水量平衡模型。
在其中一个实施例中，该待求解参数包括第一待求解参数、第二待求解参数、第三待求解参数、第四待求解参数、第五待求解参数和第六待求解参数；
该水量输出通量函数的构建过程，包括：
根据实际蒸散发量函数、土壤深层渗漏函数以及流域径流量函数联合得到该水量输出通量函数；
该实际蒸散发量函数包括流域月尺度潜在蒸散发、土壤体积含水率、该第一待求解参数以及第二待求解参数之间的函数关系；其中，第一待求解参数表示估算的该潜在蒸散发的系统偏差，第二待求解参数用于确定该实际蒸散发量函数的曲线形状；
该土壤深层渗漏函数包括土壤体积含水率、该第三待求解参数以及第四待求解参数之间的函数关系；其中，第三待求解参数表示深层渗漏量的最大值，第四待求解参数表示土壤含水层水量流失的速率；
该流域径流量函数包括该水量输入通量函数、土壤体积含水率以及该第五待求解参数之间的函数关系；其中，第五待求解参数表示所述水量输入通量函数中降水产流的速率；
该方法还包括：
根据该水量平衡模型确定第六待求解参数；其中，第六待求解参数表示该流域研究土层的厚度。
第二方面，提供了一种无资料地区径流估算装置，该装置包括：
获取模块，用于获取驱动数据，该驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；
构建模块，用于构建流域月尺度水量平衡变化的目标函数，确定该目标函数的待求解参数；
初始化模块，用于对该待求解参数进行初始化，得到待求解参数的初始值；以及将该待求解参数的初始值和该驱动数据输入该目标函数，得到目标函数的初始函数值；
优化模块，用于基于该待求解参数的初始值和该初始函数值，对该待求解参数进行迭代计算，得到待求解参数的优化值，将该待求解参数的优化值输入该目标函数得到目标函数值；
径流计算模块，用于在该目标函数值满足预设的迭代终止条件的情况下，根据该待求解参数的优化值估算预设时间段内流域月尺度径流值。
第三方面，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行该计算机程序时实现以下步骤：
获取驱动数据，该驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；
构建流域月尺度水量平衡变化的目标函数，确定该目标函数的待求解参数；
对该待求解参数进行初始化，得到待求解参数的初始值；
将该待求解参数的初始值和该驱动数据输入该目标函数，得到目标函数的初始函数值；
基于该待求解参数的初始值和该初始函数值，对该待求解参数进行迭代计算，得到待求解参数的优化值，将该待求解参数的优化值输入该目标函数得到目标函数值；
在该目标函数值满足预设的迭代终止条件的情况下，根据该待求解参数的优化值估算预设时间段内流域月尺度径流值。
第四方面，提供了一种计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现以下步骤：
获取驱动数据，该驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；
构建流域月尺度水量平衡变化的目标函数，确定该目标函数的待求解参数；
对该待求解参数进行初始化，得到待求解参数的初始值；
将该待求解参数的初始值和该驱动数据输入该目标函数，得到目标函数的初始函数值；
基于该待求解参数的初始值和该初始函数值，对该待求解参数进行迭代计算，得到待求解参数的优化值，将该待求解参数的优化值输入该目标函数得到目标函数值；
在该目标函数值满足预设的迭代终止条件的情况下，根据该待求解参数的优化值估算预设时间段内流域月尺度径流值。
上述无资料地区径流估算方法、装置、计算机设备和存储介质，通过获取驱动数据，所述驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；然后，构建流域月尺度水量平衡变化的目标函数，并确定目标函数的待求解参数；然后对待求解参数进行初始化，并将待求解参数的初始值与驱动数据输入目标函数，得到目标函数的初始函数值；并基于所述待求解参数的初始值和所述初始函数值，对所述待求解参数进行迭代计算；最后，根据所述待求解参数的优化值估算预设时间段内流域月尺度径流值。使用该方法可以在不依赖任何实测径流资料的条件下，通过输入较少的驱动数据即可对无资料地区的径流值进行估算，估算准确率高。
附图说明
图1为一个实施例中无资料地区径流估算方法的流程示意图；
图2为一个实施例中流域月尺度水量平衡模型的示意图；
图3为一个实施例中无资料地区径流估算方法的总体框架图；
图4为一个实施例中无资料地区径流估算方法的测试对象图；
图5为一个实施例中无资料地区径流估算方法的测试结果图；
图6为一个实施例中无资料地区径流估算装置的结构框图；
图7为一个实施例中计算机设备的内部结构图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
河流径流量是水文循环中的重要变量，也是用于区域水资源评估、生态环境保护和气候变化分析的基础数据。依据径流形成的水文过程，总径流量包括地表径流，壤中流和地下径流。相关的径流监测数据主要来源于水文监测站点，虽然其精度较高，但是全球站点分布十分不均匀且多数站点监测时段较短，受地形及观测成本的制约，许多无资料区缺少可靠的径流数据。水文学研究中通常采用搭建水文模型的方法模拟流域径流量，然而相关方法依赖实测径流资料对模型参数进行率定，在无资料地区应用受限。
无资料地区的径流模拟问题是水文学研究的难点之一，常用的方法主要有比拟法、插值法、回归法等参数移植法。其中，比拟法和插值法需要使用相邻流域的实测流量数据，且对模拟流域和参考流域之间的地理位置和地形空间变换很敏感；回归法通过建立参考流域的流域特征和模型参数之间的回归关系，然后将这种关系用于其他流域，但其效果往往很有限。由于不同流域的独特性，参数移植的方法具有很大的不确定性。
近年来，基于水量平衡公式，由流域的状态变量(土壤水储量)推求各通量(降水、蒸散发、径流、侧向流和下渗量、灌溉量)的方法逐渐发展，例如Brocca等利用实测土壤含水量变化估算降水量，Akbar等利用遥感卫星SMAP观测的土壤水分信息估算网格尺度的蒸散发及排泄量，Filippucci等利用实测土壤含水量及实测降水信息推算灌区的灌溉量。这类方法物理原理明晰，相比于传统水文模型而言输入的驱动数据较少，求解方法多采用数学优化思想，可以克服水文模型中依赖实测资料进行参数率定的局限。然而，目前这类研究通常在站点或网格尺度开展，假设径流量可忽略，例如Brocca等，Akbar等在研究中假设所有降水量均可入渗，不产生地表径流；Filippucci等假设在计算灌溉量时地表径流量可忽略，因为优化的灌溉方式不期望产生多余的水量流出。这些假设可应用于较小尺度和特定情境，但是在流域尺度上，忽略径流量会导致水量平衡不闭合，引起较大误差。总之，此类研究方法目前还未应用于流域径流量的估算中，其在径流过程解析、水文循环研究及水资源管理方面的潜力有待进一步发掘。
为了获取无资料地区径流序列，本申请基于流域月尺度水量平衡模型，利用高维参数空间的优化求解思想，提供一种不依赖任何实测径流资料的月尺度径流求解方法，可为无资料地区的径流监测提供参考。需要说明的是，本申请所指无资料地区是指研究流域没有实测的径流资料，但是研究流域可能会有其他少量的观测资料，例如流域月尺度降水等，以这些观测资料作为参考，可以加快本申请目标函数的收敛速度，研究流域没有观测资料也不影响最终求解的径流值。本申请通过全球再分析数据集ERA5获取驱动数据，对无资料地区的径流值进行估算。该方法成功应用于青藏高原江河源区，包括雅鲁藏布江(以奴下站为流域出水口)、怒江(以嘉玉桥站为流域出水口)、澜沧江(以昌都站为流域出水口)、长江(以直门达站为流域出水口)、黄河(以唐乃亥站为流域出水口)源区，可为全球其他无资料流域的径流估算提供借鉴。
在一个实施例中，如图1所示，提供了一种无资料地区径流估算方法，本实施例以该方法应用于终端进行举例说明，可以理解的是，该方法也可以应用于服务器，还可以应用于包括终端和服务器的系统，并通过终端和服务器的交互实现。本实施例中，该方法包括以下步骤：
步骤102，获取驱动数据，该驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率。
其中，驱动数据是用来进行流域月尺度水量平衡变化的目标函数求解的必要参数，以对流域月尺度的径流值进行估算。预设时间段可以根据求解的需求来确定，本申请实施例对此不做限定，例如要估算2019年1-12月的径流，预设时间段就是2019年1-12月，则获取对应的时间段内的驱动数据。
具体地，从全球再分析数据集ERA5中获取预设时间段内的驱动数据，本申请实施例中的驱动数据包括流域月尺度降水、潜在蒸散发和土壤体积含水率等。ERA5是第五代ECMWF(European Centre for Medium-Range Weather Forecasts，欧洲中期天气预报中心)大气再分析全球气候数据，ERA5数据提供每小时的大气、陆地和海洋气候变量的估计值，地球数据精确到了30km网格，包括了137层的大气数据。ERA5数据集为netcdf格式，具体可利用Matlab程序中“ncread”函数读取。
步骤104，构建流域月尺度水量平衡变化的目标函数，确定该目标函数的待求解参数。
具体地，根据两种计算方式得到研究流域土壤水储量变化关系，接着根据这两种方式计算的土壤水储量变化的均方根误差构建需要优化的目标函数，并确定目标函数中的待求解参数，构成参数空间。可以根据由水量输入以及水量输出通量计算的土壤水储量变化与由研究流域月尺度水量状态量的变化之间的均方根误差来构建目标函数，将目标函数中的未知参数作为待求解参数。
步骤106，对该待求解参数进行初始化，得到待求解参数的初始值。
具体地，为了求解目标函数，需要对待求解参数赋初始值，即对上述参数空间赋初值，得到参数空间的初始值。理论上初始值可以任意赋值，但是实际求解中，不同的初始值会影响目标函数的优化曲线收敛的快慢程度，为了避免无效运算，可以参考研究流域或相邻流域的实际情况，结合不同参数的物理意义进行初始值赋值，相对合理的初始值可以使求解过程更快收敛，但不影响最终的优化结果。例如，假设用a表示ERA5估算的潜在蒸散发的系统偏差值，则可以用从ERA5中获取的潜在蒸散发与观测记录中的蒸散发的差值给a赋初值。
需要明确的是，观测记录可以帮助确定参数的初始值，使研究流域目标函数的初始状态更贴近真实状态，以加快目标函数的收敛过程，不会影响最终求解的径流值。因此，本申请实施例对研究流域是否有其他变量的观测资料不做要求。
步骤108，将该待求解参数的初始值和该驱动数据输入该目标函数，得到目标函数的初始函数值。
具体地，将上述参数空间的初始值与获取的驱动数据一起输入目标函数进行计算，得到目标函数的初始函数值。初始函数值用于跟后续优化后的目标函数值进行比较，从而判断是否达到迭代终止条件。
步骤110，基于该待求解参数的初始值和该初始函数值，对该待求解参数进行迭代计算，得到待求解参数的优化值，将该待求解参数的优化值输入该目标函数得到目标函数值。
具体地，可以利用数学分析中的最优化方法对待求解参数进行迭代优化，得到待求解参数的优化值，即参数空间的优化值。例如，可以使用共轭梯度法对目标函数中的待求解参数进行优化求解，还可以用其他最优化方法进行求解，本实施例在此不做限定。最后，将该待求解参数的优化值输入该目标函数得到目标函数值。
步骤112，在该目标函数值满足预设的迭代终止条件的情况下，根据该待求解参数的优化值估算预设时间段内流域月尺度径流值。
具体地，终止条件是对目标函数求解精度的要求。在本申请中设置优化阈值为土壤水储量变化量级的1‰，终止条件为：当两次优化的目标函数值(即本轮优化得到的目标函数值与上一轮优化得到的目标函数值)的差值小于优化阈值时，迭代终止；若两次优化的目标函数值的差值大于或者等于优化阈值时，则继续下一轮的迭代计算。再根据本轮待求解参数的优化值求解预设时间段内流域月尺度径流值。
例如，土壤水储量变化量级用
Figure BDA0002901973590000091
表示，本轮优化得到的目标函数值用Ji+1表示，上一轮优化得到的目标函数值用Ji表示，当两次优化的效果提升小于优化阈值时，认为找到最优解，即当
Figure BDA0002901973590000092
时，优化过程终止。最后，根据本轮优化的待求解参数的参数空间Xi+1以及径流函数求解流域月尺度径流值。式中
Figure BDA0002901973590000101
表示参与求解的所有样本的水储量变化的绝对值的均值。
上述无资料地区径流估算方法中，通过获取驱动数据，所述驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；然后，构建流域月尺度水量平衡变化的目标函数，并确定目标函数的待求解参数；然后对待求解参数进行初始化，并将待求解参数的初始值与驱动数据输入目标函数，得到目标函数的初始函数值；并基于所述待求解参数的初始值和所述初始函数值，对所述待求解参数进行迭代计算；最后，根据所述待求解参数的优化值估算预设时间段内流域月尺度径流值。使用该方法无需实测的径流数据即可对无资料地区的径流值进行估算，估算准确率高。
在一个实施例中，对该待求解参数进行迭代计算，得到待求解参数的优化值，包括：
计算该目标函数对该待求解参数的梯度函数；
计算该梯度函数在上一轮待求解参数的优化值下的梯度值；
计算本轮该待求解参数在该梯度值下的参数值，得到该待求解参数的优化值。
具体地，梯度函数表示目标函数对待求解参数的梯度，可以将目标函数对待求解参数进行求导，例如，可利用Matlab程序中“diff”函数求得目标函数对待求解参数的梯度。例如，目标函数用J表示，本轮待求解参数用xi+1表示，上一轮待求解参数用xi表示；目标函数对待求解参数的梯度则用
Figure BDA0002901973590000102
来表示。
然后求目标函数对待求解参数的梯度在上一轮待求解参数的优化值下的梯度值。具体地，将上一轮待求解参数xi带入梯度函数中，即可得到待求解参数的梯度在xi条件下的梯度值。
最后，根据上一轮待求解参数xi、上述梯度值进行优化计算，得到本轮待求解参数的优化值xi+1，根据本轮待求解参数的优化值输入目标函数计算，得到本轮的目标函数值Ji+1。
在一个实施例中，计算本轮该待求解参数在该梯度值下的参数值，得到该待求解参数的优化值，包括：
根据该上一轮待求解参数的优化值、该梯度值以及预设的优化精度，将所述上一轮待求解参数的优化值沿着所述梯度值的负方向进行迭代计算，得到该待求解参数的优化值。
具体地，根据该上一轮待求解参数的优化值、该梯度值以及预设的优化精度将待求解参数沿梯度负方向进行移动。在有最优解(有极值)的高维凸空间中，参数空间梯度的反方向是靠近极值的最快的方向，因此每次沿着参数空间的反方向进行移动，如
Figure BDA0002901973590000111
就能逼近最优解，即待求解参数的优化值，此处α代表优化精度，优化精度可以任意给定，本申请实施例中设置为5，也可以根据实际需求设置为其他数值，需要根据运算和求解的要求给出合理值。α越小，精度越高，但是计算成本越高，实际操作中α一般可取1-10。
在一个实施例中，该待求解参数包括多个待求解参数，共同组成参数空间；
计算该目标函数对该待求解参数的梯度函数，包括：
计算该目标函数对该参数空间中每一个待求解参数的梯度函数；
计算该梯度函数在上一轮待求解参数的优化值下的梯度值，包括：
计算该每一个待求解参数的梯度函数在对应的上一轮待求解参数的优化值下的梯度值；
根据该上一轮待求解参数的优化值、该梯度值以及优化精度，将该上一轮待求解参数的优化值沿着该梯度值的负方向进行迭代计算，得到该待求解参数的优化值，包括：
根据该上一轮每一个待求解参数的优化值、对应的该梯度值以及优化精度将该上一轮每一个待求解参数的优化值沿着该梯度值的负方向进行迭代计算，得到对应的优化后的参数值，该每一个优化后的参数值共同组成优化的参数空间。
具体地，目标函数中包含多个待求解参数，共同组成参数空间，例如，可以用X表示待求参数的参数空间，X0表示参数空间的初始值，在迭代优化的过程中需要对每一个待求解参数分别进行优化计算，进而得到每一个待求解参数的优化值。首先计算目标函数对参数空间中每一个待求解参数的梯度函数，然后，计算每一个待求解参数的梯度函数在对应的上一轮待求解参数的优化值下的梯度值，最后根据上一轮每一个待求解参数的优化值、对应的梯度值以及优化精度将待求解参数沿梯度负方向进行移动，得到对应的优化后的参数值，该每一个优化后的参数值共同组成本轮优化的参数空间Xi+1，Xi+1＝[ai+1,bi+1,ci+1,di+1,fi+1,zi+1]。
在一个实施例中，水量平衡模型的构建过程包括：
根据预设时间段内水量输入通量函数、水量输出通量函数以及水量的状态量之间的关系构建该水量平衡模型。
具体地，根据流域月尺度水量输入通量、输出通量及状态量变化三者之间的关系构建流域月尺度水量平衡模型。水量平衡公式可表达为：一段时间内研究流域内土壤水储量的变化为水量输入(降水)与输出(实际蒸散发、径流和深层渗漏)的差值。
如图2所示，以一定厚度土层为研究对象，输入通量为流域降水量，输出通量包括流域实际蒸散发、径流(包括地表径流、壤中流及地下径流)及土壤深层渗漏，流域状态量的变化为该厚度的土壤水储量变化。例如，流域月尺度的水量平衡关系式可表示为：
Figure BDA0002901973590000121
式中，Δz表示研究流域土层的厚度，在本申请中作为待求参数之一，是满足水量平衡闭合的参数，Δz不仅包括实测土壤水分的土层厚度，也与更深层的土壤水分动态有关；θ表示土壤的体积含水率；P(t)表示月降水量，单位为mm/mon，即流域的水量输入通量；L(θ)表示水量输出通量，单位为mm/mon。
本申请实施例中流域月尺度下的θ值和P(t)值可由全球再分析数据集ERA5获取。Δz是自由参数，需要赋一个初始值，然后在参数迭代的过程中不断更新，最终找到最优解。实际操作中，为了让参数迭代收敛更快，给定初始值时，可以参考ERA5的土壤层厚度(289cm)给定一个约300cm的初始值。
在一个实施例中，构建流域月尺度水量平衡变化的目标函数，包括：
根据水量平衡模型中水量的通量变化关系进行转换，得到第一土壤水储量变化参数；
根据该水量平衡模型中水量的状态量变化关系，得到第二土壤水储量变化参数；
根据该第一土壤水储量变化参数、该第二土壤水储量变化参数构建该流域月尺度水量平衡变化的目标函数。
具体地，根据水量平衡模型中水量的输入通量、输出通量的差值以及研究流域土层的厚度计算得到第一土壤水储量变化参数。根据研究流域水量的状态量计算的土壤水储量变化作为第二土壤水储量变化参数。最后根据第一土壤水储量变化参数、该第二土壤水储量变化参数的均方根误差构建该流域月尺度水量平衡变化的目标函数。如下列公式(2)所示：
Figure BDA0002901973590000131
式中N表示样本数目，X表示待求参数的参数空间，如有a,b,c,d,f,z六个待求解参数，则X＝[a,b,c,d,f,z]。
Figure BDA0002901973590000132
表示第一土壤水储量变化参数，
Figure BDA0002901973590000133
表示第二土壤水储量变化参数。样本数目N是参与求解的数据的个数，根据待求解的时段确定，例如估算2019年1-12月的径流值，样本数目为12。
在一个实施例中，该待求解参数包括第一待求解参数、第二待求解参数、第三待求解参数、第四待求解参数、第五待求解参数和第六待求解参数；
该水量输出通量函数的构建过程，包括：
根据实际蒸散发量函数、土壤深层渗漏函数以及流域径流量函数联合得到该水量输出通量函数；
该实际蒸散发量函数包括流域月尺度潜在蒸散发、土壤体积含水率、该第一待求解参数以及第二待求解参数之间的函数关系；其中，第一待求解参数表示估算的该潜在蒸散发的系统偏差，第二待求解参数用于确定该实际蒸散发量函数的曲线形状；
该土壤深层渗漏函数包括土壤体积含水率、该第三待求解参数以及第四待求解参数之间的函数关系；其中，第三待求解参数表示深层渗漏量的最大值，第四待求解参数表示土壤含水层水量流失的速率；
该流域径流量函数包括该水量输入通量函数、土壤体积含水率以及该第五待求解参数之间的函数关系；其中，第五待求解参数表示所述水量输入通量函数中降水产流的速率；
该方法还包括：
根据该水量平衡模型确定第六待求解参数；其中，第六待求解参数表示该流域研究土层的厚度。
具体地，根据实际蒸散发量函数、土壤深层渗漏函数以及流域径流量函数联合得到该水量输出通量函数。例如，水量输出通量函数L(θ)包含的分量如下式所示：
L(θ)＝ET(θ,a,b)+D(θ,c,d)+R(θ,f) (3)
式中，ET(θ,a,b)表示实际蒸散发，单位是mm/mon；D(θ,c,d)表示土壤深层渗漏，单位为mm/mon；R(θ,f)表示流域径流深，由地表径流和地下径流组成，单位为mm/mon。a,b,c,d,f均为确定各分量的数学函数形式的待求参数，分别为第一待求解参数、第二待求解参数、第三待求解参数、第四待求解参数、第五待求解参数。
公式(3)中的各分量均可用数学函数表达，其中实际蒸散发ET可由双曲正切函数拟合(公式(4))，其渐近线为潜在蒸散发的量值(表征流域蒸发能力的最大值)，关键参数需确定函数曲线形状，表征实际蒸散发受水分胁迫向受能量胁迫的转化过程。该函数的物理意义为：当土壤体积含水率为零时，流域的实际蒸散发为零，此时实际蒸散发受水分胁迫；当土壤体积含水率较大时，即土壤足够湿润时，流域的实际蒸散发接近潜在蒸散发，函数值逼近渐近线值域，此时实际蒸散发受能量胁迫。
Figure BDA0002901973590000151
式中，PET(t)表示流域月尺度潜在蒸散发，单位为mm/mon，可由全球再分析数据集ERA5获取；第一待求解参数a表示ERA5估算的潜在蒸散发的系统偏差；关于第二待求解参数b的函数sig(b)为sigmoid函数，即
Figure BDA0002901973590000152
Figure BDA0002901973590000153
表示土壤孔隙度，可由土壤中砂土和黏土的百分比计算，即
Figure BDA0002901973590000154
Figure BDA0002901973590000155
sand，clay分别表示砂土、黏土在总土层中的质量百分比，该数据可由全球土壤性质数据集获取。
在初始化阶段，将土壤孔隙度
Figure BDA0002901973590000158
与流域月尺度降水、潜在蒸散发和土壤体积含水率等驱动数据以及待求解参数的初始值一起输入目标函数，得到目标函数的初始函数值。
深层渗漏量D包括土壤侧向流和垂向深层渗漏两部分，其数学函数形式可表达为土壤体积含水率的幂函数：
Figure BDA0002901973590000156
公式(5)中第三待求解参数c可表示深层渗漏量的最大值，第四待求解参数d表示土壤含水层水量流失的速率。
径流R的数学函数形式可表达为考虑降水和土壤体积含水率的幂函数，其最大值等于流域降水量：
Figure BDA0002901973590000157
公式(6)中第五待求解参数f可表示降水产流的速率，即降水转化为径流的速率。
水量平衡模型(即公式(1))中第六待求解参数Δz表示研究流域土层的厚度。
与传统水文模型依赖大量地面实测径流量率定模型参数的方法不同，本申请实施例所提供的无资料地区径流估计方法不需要任何实测径流资料，所需驱动数据少，公式概化简洁，适用性强，可服务于水资源利用与管理、水文过程及气候变化研究等诸多领域。使用该方法可为无资料地区径流监测提供数据参考，有效解决水文模型对实测径流资料的要求。
本申请实施例还将所提供无资料地区径流估算方法在青藏高原五个江河源区进行测试，模拟了2000-2017年流域月尺度径流量，利用雅鲁藏布江(以奴下站为流域出水口)、怒江(以嘉玉桥站为流域出水口)、澜沧江(以昌都站为流域出水口)、长江(以直门达站为流域出水口)、黄河(以唐乃亥站为流域出水口)源区流域出水口水文站观测的径流量进行检验。本申请实施例的检验指标为纳什效率系数NSE，计算公式如下：
Figure BDA0002901973590000161
式中，Qm和Qo分别是求解的径流量和观测的径流量，
Figure BDA0002901973590000162
表示观测值的平均值，该指标的理论最优值为1。
需要说明的是，虽然上述五个流域有实测的径流数据，但是本申请在测试的过程中并未利用实测的径流数据参与计算，而是将这些流域视为无资料地区，通过全球再分析数据集ERA5获取驱动数据，对这些无资料地区的径流进行估计。最后用这五个流域的实测径流数据与通过本申请估算的径流数据对本申请所提供的无资料地区径流估算方法的准确率进行检验，以验证方法的有效性。
研究流域及观测站的地理位置如图4所示，结果表明(如图5所示)：在完全不依赖实测径流的情况下，五个测试流域的月尺度径流求解结果与实测径流相比纳什效率系数(NSE)均可大于0.6，NSE在怒江流域最高(0.85)，利用优化方法求解的径流可较为准确地反映年际变化、峰值及量级。验证结果表明：本申请实施例所提供的无资料地区径流估计方法可以解决无资料地区径流模拟难题，在无资料地区适应性较强，且准确性较高。
为了易于理解本申请实施例提供的技术方案，如图3所示，以完整的无资料地区径流估计过程对本申请实施例提供的无资料地区径流估计方法进行简要说明：
(1)根据预设时间段内水量输入通量函数P(t)与水量输出通量函数L(θ)的差值以及水量的状态量
Figure BDA0002901973590000171
之间的关系构建所述水量平衡模型，即公式(1)。
(2)根据实际蒸散发量函数ET(θ,a,b)、土壤深层渗漏函数D(θ,c,d)以及流域径流量函数R(θ,f)联合得到所述水量输出通量函数，即公式(3)。
所述实际蒸散发量函数(即公式(4))包括流域月尺度潜在蒸散发、土壤体积含水率、所述第一待求解参数以及第二待求解参数之间的函数关系；其中，第一待求解参数表示估算的所述潜在蒸散发的系统偏差，第二待求解参数用于确定所述实际蒸散发量函数的曲线形状；
所述土壤深层渗漏函数(即公式(5))包括土壤体积含水率、所述第三待求解参数以及第四待求解参数之间的函数关系；其中，第三待求解参数表示深层渗漏量的最大值，第四待求解参数表示土壤含水层水量流失的速率；
所述流域径流量函数(即公式(6))包括所述水量输入通量函数、土壤体积含水率以及所述第五待求解参数之间的函数关系；其中，第五待求解参数表示所述水量输入通量函数中降水产流的速率；
根据所述水量平衡模型确定第六待求解参数；其中，第六待求解参数表示所述流域研究土层的厚度。
(3)根据水量平衡模型中水量的通量变化关系进行转换，得到第一土壤水储量变化参数；
根据所述水量平衡模型中水量的状态量变化关系，得到第二土壤水储量变化参数；
根据所述第一土壤水储量变化参数、所述第二土壤水储量变化参数的均方根误差构建所述流域月尺度水量平衡变化的目标函数，即公式(2)。
(4)获取驱动数据，所述驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率。
(5)对所述待求解参数进行初始化，得到待求解参数的初始值；
将所述待求解参数的初始值和所述驱动数据输入所述目标函数，得到目标函数的初始函数值。
(6)基于所述待求解参数的初始值和所述初始函数值，对所述待求解参数进行迭代计算，即计算所述目标函数对所述多个待求解参数中每一个待求解参数的梯度函数，计算所述每一个待求解参数的梯度函数在对应的上一轮待求解参数的优化值下的梯度值，根据所述上一轮每一个待求解参数的优化值、对应的所述梯度值以及优化精度进行计算，得到对应的优化后的参数值，所述每一个优化后的参数值共同组成优化的参数空间，将所述待求解参数的优化值输入所述目标函数得到目标函数值。
(7)在所述目标函数值满足预设的迭代终止条件的情况下，根据所述待求解参数的优化值估算流域月尺度径流值。
(8)对估算的流域月尺度径流值的准确率进行计算。
应该理解的是，虽然图1、3的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1、3中的至少一部分步骤可以包括多个步骤或者多个阶段，这些步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤中的步骤或者阶段的至少一部分轮流或者交替地执行。
在一个实施例中，如图6所示，提供了一种无资料地区径流估算装置，包括：获取模块602、构建模块604、初始化模块606、优化模块608和径流计算模块610，其中：
获取模块602，用于获取驱动数据，该驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率。
构建模块604，用于构建流域月尺度水量平衡变化的目标函数，确定该目标函数的待求解参数。
初始化模块606，用于对该待求解参数进行初始化，得到待求解参数的初始值；以及将该待求解参数的初始值和该驱动数据输入该目标函数，得到目标函数的初始函数值。
优化模块608，用于基于该待求解参数的初始值和该初始函数值，对该待求解参数进行迭代计算，得到待求解参数的优化值，将该待求解参数的优化值输入该目标函数得到目标函数值。
径流计算模块610，用于在该目标函数值满足预设的迭代终止条件的情况下，根据该待求解参数的优化值估算预设时间段内流域月尺度径流值。
在一个实施例中，构建模块604还用于根据预设时间段内水量输入通量函数、水量输出通量函数以及水量的状态量之间的关系构建该水量平衡模型。
在一个实施例中，构建模块604还用于根据水量平衡模型中水量的通量变化关系进行转换，得到第一土壤水储量变化参数；
以及根据该水量平衡模型中水量的状态量变化关系，得到第二土壤水储量变化参数；
以及根据该第一土壤水储量变化参数、该第二土壤水储量变化参数构建该流域月尺度水量平衡变化的目标函数。
在一个实施例中，该待求解参数包括第一待求解参数、第二待求解参数、第三待求解参数、第四待求解参数、第五待求解参数和第六待求解参数；
构建模块604还用于根据实际蒸散发量函数、土壤深层渗漏函数以及流域径流量函数联合得到该水量输出通量函数；
该实际蒸散发量函数包括流域月尺度潜在蒸散发、土壤体积含水率、该第一待求解参数以及第二待求解参数之间的函数关系；其中，第一待求解参数表示估算的该潜在蒸散发的系统偏差，第二待求解参数用于确定该实际蒸散发量函数的曲线形状；
该土壤深层渗漏函数包括土壤体积含水率、该第三待求解参数以及第四待求解参数之间的函数关系；其中，第三待求解参数表示深层渗漏量的最大值，第四待求解参数表示土壤含水层水量流失的速率；
该流域径流量函数包括该水量输入通量函数、土壤体积含水率以及该第五待求解参数之间的函数关系；其中，第五待求解参数表示所述水量输入通量函数中降水产流的速率；
构建模块604还用于根据该水量平衡模型确定第六待求解参数；其中，第六待求解参数表示该流域研究土层的厚度。
在一个实施例中，无资料地区径流估算装置还包括验证模块，验证模块用于对估算的流域月尺度径流值的准确率进行计算。
在一个实施例中，优化模块608还用于计算该目标函数对该待求解参数的梯度函数；
以及计算该梯度函数在上一轮待求解参数的优化值下的梯度值；
以及计算本轮该待求解参数在该梯度值下的参数值，得到该待求解参数的优化值。
在一个实施例中，优化模块608还用于根据该上一轮待求解参数的优化值、该梯度值以及预设的优化精度，将所述上一轮待求解参数的优化值沿着所述梯度值的负方向进行迭代计算，得到该待求解参数的优化值。
在一个实施例中，该待求解参数包括多个待求解参数，共同组成参数空间；
优化模块608还用于计算该目标函数对该参数空间中每一个待求解参数的梯度函数；
以及计算该每一个待求解参数的梯度函数在对应的上一轮待求解参数的优化值下的梯度值；
以及根据该上一轮每一个待求解参数的优化值、对应的该梯度值以及优化精度，将该上一轮每一个待求解参数的优化值沿着该梯度值的负方向进行迭代计算，得到对应的优化后的参数值，该每一个优化后的参数值共同组成优化的参数空间。
关于无资料地区径流估算装置的具体限定可以参见上文中对于无资料地区径流估算方法的限定，在此不再赘述。上述无资料地区径流估算装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图7所示。该计算机设备包括通过系统总线连接的处理器、存储器、通信接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的通信接口用于与外部的终端进行有线或无线方式的通信，无线方式可通过WIFI、运营商网络、NFC(近场通信)或其他技术实现。该计算机程序被处理器执行时以实现一种无资料地区径流估算方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图7中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，还提供了一种计算机设备，包括存储器和处理器，存储器中存储有计算机程序，该处理器执行计算机程序时实现上述各方法实施例中的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现上述各方法实施例中的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和易失性存储器中的至少一种。非易失性存储器可包括只读存储器(Read-Only Memory，ROM)、磁带、软盘、闪存或光存储器等。易失性存储器可包括随机存取存储器(Random Access Memory，RAM)或外部高速缓冲存储器。作为说明而非局限，RAM可以是多种形式，比如静态随机存取存储器(Static Random Access Memory，SRAM)或动态随机存取存储器(Dynamic Random Access Memory，DRAM)等。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种无资料地区径流估算方法，其特征在于，所述方法包括：
获取驱动数据，所述驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；
根据由水量输入以及水量输出通量计算的土壤水储量变化与由研究流域月尺度水量状态量的变化之间的均方根误差来构建流域月尺度水量平衡变化的目标函数，将所述目标函数中的未知参数作为待求解参数；
对所述待求解参数进行初始化，得到待求解参数的初始值；
将所述待求解参数的初始值和所述驱动数据输入所述目标函数，得到目标函数的初始函数值；
基于所述待求解参数的初始值和所述初始函数值，计算所述目标函数对所述待求解参数的梯度函数；计算所述梯度函数在上一轮待求解参数的优化值下的梯度值；根据所述上一轮待求解参数的优化值、所述梯度值以及预设的优化精度，将所述上一轮待求解参数的优化值沿着所述梯度值的负方向进行迭代计算，得到所述待求解参数的优化值；将所述待求解参数的优化值输入所述目标函数得到目标函数值；
在所述目标函数值满足预设的迭代终止条件的情况下，根据所述待求解参数的优化值估算预设时间段内流域月尺度径流值。
2.根据权利要求1所述的方法，其特征在于，所述待求解参数包括多个待求解参数，共同组成参数空间；
所述计算所述目标函数对所述待求解参数的梯度函数，包括：
计算所述目标函数对所述参数空间中每一个待求解参数的梯度函数；
所述计算所述梯度函数在上一轮待求解参数的优化值下的梯度值，包括：
计算所述每一个待求解参数的梯度函数在对应的上一轮待求解参数的优化值下的梯度值；
所述根据所述上一轮待求解参数的优化值、所述梯度值以及优化精度，将所述上一轮待求解参数的优化值沿着所述梯度值的负方向进行迭代计算，得到所述待求解参数的优化值，包括：
根据所述上一轮每一个待求解参数的优化值、对应的所述梯度值以及优化精度，将所述上一轮每一个待求解参数的优化值沿着所述梯度值的负方向进行迭代计算，得到对应的优化后的参数值，所述每一个优化后的参数值共同组成优化的参数空间；
所述将所述待求解参数的优化值输入所述目标函数得到目标函数值，包括：
将所述优化的参数空间输入所述目标函数得到目标函数值。
3.根据权利要求1所述的方法，其特征在于，所述构建流域月尺度水量平衡变化的目标函数，包括：
根据水量平衡模型中水量的输入通量与输出通量的差值除以研究流域土层的厚度得到第一土壤水储量变化参数；
根据研究流域水量的状态量计算得到的土壤水储量变化作为第二土壤水储量变化参数；
根据所述第一土壤水储量变化参数、所述第二土壤水储量变化参数的均方根误差构建所述流域月尺度水量平衡变化的目标函数。
4.根据权利要求3所述的方法，其特征在于，所述水量平衡模型的构建过程包括：
根据预设时间段内水量输入通量函数、水量输出通量函数以及水量的状态量之间的关系构建所述水量平衡模型。
5.根据权利要求4所述的方法，其特征在于，所述待求解参数包括第一待求解参数、第二待求解参数、第三待求解参数、第四待求解参数、第五待求解参数和第六待求解参数；
所述水量输出通量函数的构建过程，包括：
根据实际蒸散发量函数、土壤深层渗漏函数以及流域径流量函数联合得到所述水量输出通量函数；
所述实际蒸散发量函数包括流域月尺度潜在蒸散发、土壤体积含水率、所述第一待求解参数以及第二待求解参数之间的函数关系；其中，第一待求解参数表示估算的所述潜在蒸散发的系统偏差，第二待求解参数用于确定所述实际蒸散发量函数的曲线形状；
所述土壤深层渗漏函数包括土壤体积含水率、所述第三待求解参数以及第四待求解参数之间的函数关系；其中，第三待求解参数表示深层渗漏量的最大值，第四待求解参数表示土壤含水层水量流失的速率；
所述流域径流量函数包括所述水量输入通量函数、土壤体积含水率以及所述第五待求解参数之间的函数关系；其中，第五待求解参数表示所述水量输入通量函数中降水产流的速率；
所述方法还包括：
根据所述水量平衡模型确定第六待求解参数；其中，第六待求解参数表示所述流域研究土层的厚度。
6.一种无资料地区径流估算装置，其特征在于，所述装置包括：
获取模块，用于获取驱动数据，所述驱动数据包括预设时间段内流域月尺度降水量、潜在蒸散发量和土壤体积含水率；
构建模块，用于根据由水量输入以及水量输出通量计算的土壤水储量变化与由研究流域月尺度水量状态量的变化之间的均方根误差来构建目标函数，将目标函数中的未知参数作为待求解参数；
初始化模块，用于对所述待求解参数进行初始化，得到待求解参数的初始值；以及将所述待求解参数的初始值和所述驱动数据输入所述目标函数，得到目标函数的初始函数值；
优化模块，用于基于所述待求解参数的初始值和所述初始函数值，计算所述目标函数对所述待求解参数的梯度函数；计算所述梯度函数在上一轮待求解参数的优化值下的梯度值；根据所述上一轮待求解参数的优化值、所述梯度值以及预设的优化精度，将所述上一轮待求解参数的优化值沿着所述梯度值的负方向进行迭代计算，得到所述待求解参数的优化值，将所述待求解参数的优化值输入所述目标函数得到目标函数值；
径流计算模块，用于在所述目标函数值满足预设的迭代终止条件的情况下，根据所述待求解参数的优化值估算预设时间段内流域月尺度径流值。
7.根据权利要求6所述的装置，其特征在于，所述待求解参数包括多个待求解参数，共同组成参数空间；
优化模块还用于计算所述目标函数对所述参数空间中每一个待求解参数 的梯度函数；
以及计算所述每一个待求解参数的梯度函数在对应的上一轮待求解参数的优化值下的梯度值；
以及根据所述上一轮每一个待求解参数的优化值、对应的该梯度值以及优化精度，将所述上一轮每一个待求解参数的优化值沿着该梯度值的负方向进行迭代计算，得到对应的优化后的参数值，所述每一个优化后的参数值共同组成优化的参数空间。
8.根据权利要求6所述的装置，其特征在于，所述构建模块还用于根据水量平衡模型中水量的输入通量与输出通量的差值除以研究流域土层的厚度得到第一土壤水储量变化参数；
以及根据研究流域水量的状态量计算得到的土壤水储量变化作为第二土壤水储量变化参数；
以及根据所述第一土壤水储量变化参数、所述第二土壤水储量变化参数的均方根误差构建所述流域月尺度水量平衡变化的目标函数。
9.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至5中任一项所述的方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至5中任一项所述的方法的步骤。
CN202110065530.4A 2021-01-18 2021-01-18 无资料地区径流估算方法、装置、计算机设备和存储介质 Active CN113255093B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202110065530.4A CN113255093B (zh) 2021-01-18 2021-01-18 无资料地区径流估算方法、装置、计算机设备和存储介质
US17/697,915 US20230030888A1 (en) 2021-01-18 2022-03-17 Runoff estimating method and device for ungauged region, computer device, and storage medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110065530.4A CN113255093B (zh) 2021-01-18 2021-01-18 无资料地区径流估算方法、装置、计算机设备和存储介质

Publications (2)

Publication Number Publication Date
CN113255093A CN113255093A (zh) 2021-08-13
CN113255093B true CN113255093B (zh) 2022-02-15

Family

ID=77180779

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110065530.4A Active CN113255093B (zh) 2021-01-18 2021-01-18 无资料地区径流估算方法、装置、计算机设备和存储介质

Country Status (2)

Country Link
US (1) US20230030888A1 (zh)
CN (1) CN113255093B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239274A (zh) * 2021-12-17 2022-03-25 水利部交通运输部国家能源局南京水利科学研究院 多源遥感数据驱动指数滤波模型推求根层土壤水的方法
CN116415755B (zh) * 2023-05-26 2023-11-21 江西省水利科学院(江西省大坝安全管理中心、江西省水资源管理中心) 一种基于条件概率分布的径流变异程度评估方法
CN116932990B (zh) * 2023-09-14 2024-01-09 南方科技大学 一种基于转换函数方法的水文地质参数估计方法及系统
CN117077543B (zh) * 2023-10-13 2024-01-09 中国水利水电科学研究院 基于梯度搜索洗牌复形进化算法的水文模型参数率定方法
CN117057174B (zh) * 2023-10-13 2024-01-26 长江三峡集团实业发展(北京)有限公司 一种缺资料地区径流预测的方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106971034A (zh) * 2017-03-22 2017-07-21 中国水利水电科学研究院 一种基于集合卡尔曼滤波的对无资料地区径流量推求方法
CN109035105A (zh) * 2018-06-15 2018-12-18 河海大学 一种月尺度蒸散发量的定量估算方法
CN109815305A (zh) * 2019-01-18 2019-05-28 三峡大学 一种无资料地区场次洪水径流过程反演的方法
CN110570021A (zh) * 2019-08-15 2019-12-13 清华大学 径流模拟方法、装置以及计算机设备
CN110991702A (zh) * 2019-11-13 2020-04-10 清华大学 山区降水量的计算方法、装置、计算机设备和存储介质
CN112052559A (zh) * 2020-07-30 2020-12-08 浙江大学 一种基于月水量平衡的三参数水文模型及其构建方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108614948A (zh) * 2018-05-12 2018-10-02 西北农林科技大学 大尺度区域实际蒸散发量的估算方法
CN111984700B (zh) * 2020-08-17 2023-04-28 中国电建集团成都勘测设计研究院有限公司 一种基于水热耦合平衡原理的月径流预测方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106971034A (zh) * 2017-03-22 2017-07-21 中国水利水电科学研究院 一种基于集合卡尔曼滤波的对无资料地区径流量推求方法
CN109035105A (zh) * 2018-06-15 2018-12-18 河海大学 一种月尺度蒸散发量的定量估算方法
CN109815305A (zh) * 2019-01-18 2019-05-28 三峡大学 一种无资料地区场次洪水径流过程反演的方法
CN110570021A (zh) * 2019-08-15 2019-12-13 清华大学 径流模拟方法、装置以及计算机设备
CN110991702A (zh) * 2019-11-13 2020-04-10 清华大学 山区降水量的计算方法、装置、计算机设备和存储介质
CN112052559A (zh) * 2020-07-30 2020-12-08 浙江大学 一种基于月水量平衡的三参数水文模型及其构建方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
An improvement in accuracy and spatiotemporal continuity of the MODIS precipitable water vapor product based on a data fusion approach;Xueying Li等;《Remote Sensing of Environment》;20200702;全文 *
基于SCS模式的月水文模型研究;王渺林等;《重庆交通学院学报》;20060430;第25卷(第2期);全文 *
新疆无资料地区泉水沟径流分析计算方法;吴红燕;《水资源与水工程学报》;20131031;第24卷(第5期);全文 *

Also Published As

Publication number Publication date
US20230030888A1 (en) 2023-02-02
CN113255093A (zh) 2021-08-13

Similar Documents

Publication Publication Date Title
CN113255093B (zh) 无资料地区径流估算方法、装置、计算机设备和存储介质
CN111651885B (zh) 一种智慧型海绵城市洪涝预报方法
CN110570021B (zh) 径流模拟方法、装置以及计算机设备
CN114997541B (zh) 一种基于数字孪生技术的城市内涝预测方法及预警平台
CN102682335B (zh) 精确确定区域对流层延迟的神经网络方法
CN112800636B (zh) 一种估算无资料地区流域地表水资源量的方法及系统
CN110619432A (zh) 一种基于深度学习的特征提取水文预报的方法
CN109559098B (zh) 一种海绵城市试点区低影响开发设施模拟方法
CN110597873A (zh) 降水数据估计方法、装置、设备及存储介质
CN114491967B (zh) 陆地水储量预测方法、装置、设备及存储介质
Pichuka et al. Development of a method to identify change in the pattern of extreme streamflow events in future climate: Application on the Bhadra reservoir inflow in India
Li et al. A new distributed karst-tunnel hydrological model and tunnel hydrological effect simulations
CN108733952B (zh) 一种基于序贯模拟的土壤含水量空间变异性三维表征方法
CN113591288A (zh) 一种基于克里格插值的土壤湿度数据预测方法及装置
CN108154193A (zh) 一种长时间序列降水数据降尺度方法
CN110852474B (zh) 一种基于决策树算法的陆地水储量预测方法、装置及设备
Lu et al. Estimation of the interaction between groundwater and surface water based on flow routing using an improved nonlinear muskingum-cunge method
Sabzevari Runoff prediction in ungauged catchments using the gamma dimensionless time-area method
CN103870699A (zh) 基于双层异步迭代策略的水动力学洪水演进模拟方法
CN115983158A (zh) 一种地下水模型与二维水动力模型松散耦合的方法
Liu et al. Effect of irrigation methods on groundwater recharge in alluvial fan area
Al-Mussawi Kriging of groundwater level-a case study of Dibdiba Aquifer in area of Karballa-Najaf
CN110837913B (zh) 一种基于极端梯度算法的陆地水储量预测方法及设备
Xie et al. Global soil moisture storage capacity at 0.5° resolution for geoscientific modelling
CN113640228A (zh) 一种基于多光谱影像的土壤湿度预测方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant