【Hydro】水文模型比较框架MARRMoT - 包含47个概念水文模型的Matlab代码

说明

MARRMoT是一个新的水文模型比较框架，允许不同概念水文模型结构之间的客观比较。该框架为47个独特的模型结构提供了Matlab代码，所有模型结构的标准化参数范围以及每个模型的强大数值实现。该框架提供了大量的文档，用户手册和几个工作流脚本，给予如何使用该框架的例子。包括以下模型：

• FLEX-Topo
• IHACRES
• GR4J
• TOPMODEL
• SIMHYD
• VIC
• LASCAM
• TCM
• TANK
• XINANJIANG
• HYMOD
• SACRAMENTO
• MODHYDROLOG
• HBV-96
• MCRM
• SMAR
• NAM
• HYCYMODEL
• GSM-SOCONT
• ECHO
• PRMS
• CLASSIC
• IHM19

MARRMoT框架中的模型期望接近它们所基于的原始模型，但差异是不可避免的。在许多情况下，一个模型存在多个版本，但这些版本并不总是很容易区分的。存在一定程度的模型名称相等性，其中单个名称用来指代同一基模型的各种不同版本。一个很好的例子是TOPMODEL，其中许多变体存在于基于地形指数的最初概念之上。MARRMoT模型往往基于任何给定模型的旧出版物，而不是新出版物（以接近原作者的“预期”模型），但我们的选择是务实的，以实现MARRMoT中可用通量和模型结构的更大多样性。每个模型的描述列出了构成该模型MARRMoT版本基础的论文。

MARRMoT框架以模块化的方式设置，并以单个的通量文件作为基本的构建块。模型文件通过指定每个模型的内部工作方式来为其定义一类对象，而超类文件则定义了所有模型通用的所有过程。图1显示了工具箱结构的示意图概述。

源代码

降雨径流模型的模块化评估

运行实例

workflow_example_1.m

这个示例工作流程包含一个单一模型对单一集水区的应用示例。它包括5个步骤：

1. 数据准备
2. 模型选择和设置(设置参数)
3. 模型求解器设置
4. 模型生成和设置
5. 模型运行
6. 输出可视化
   

workflow_example_2.m

这个示例工作流程包含一个单一模型对单一集水区的应用示例。它包括5个步骤：

1. 数据准备
2. 模型选择和设置(默认参数)
3. 模型求解器设置
4. 模型生成和设置
5. 模型运行
6. 输出可视化
   

workflow_example_3.m

这个示例工作流程包含多个模型对单一集水区的应用示例。它包括6个步骤：

1. 数据准备
2. 模型选择和设置
3. 模型求解器设置
   对于列表中的每个模型：
   5. 模型生成和设置
   6. 模型运行
   7. 输出可视化
   

workflow_example_4.m

这个示例工作流程包含了一个单一模型对单一集水区的校准应用示例。它包括7个步骤：

1. 数据准备
2. 模型选择和设置
3. 模型求解器设置和时间步进方案
4. 校准设置
5. 模型校准
6. 校准结果评估
7. 输出可视化

测试

原始代码根据作者描述，在MATLAB版本9.11.0.1873467（R2021b）上开发的，并使用Octave 6.4.0进行了测试。笔者使用R2017a时，出现了几个可能是兼容性的问题。

1、 结构体兼容性问题

动态字段或属性名称必须为字符矢量。
将将string转换为char，即可使用结构体动态赋值。

2、append的兼容性问题

错误使用 append (line 38)
Wrong number of input arguments for obsolete matrix-based syntax.

在MATLAB中，append函数是用来将矩阵或数组的元素添加到另一个矩阵或数组的末尾。如果你想将字符串连接到另一个字符串，可以使用cat函数或者简单的加号+。

3、修改后的MARRMoT_model.m

classdef MARRMoT_model < handle
% Superclass for all MARRMoT models% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.properties% attribute to store whether we are running MATLAB or OctaveisOctave          % 1 if we're on Octave, 0 if MATLAB% static attributes, set for each models in the model definitionnumStores         % number of model storesnumFluxes         % number of model fluxesnumParams         % number of model parametersparRanges         % default parameter rangesJacobPattern      % pattern of the Jacobian matrix of model store ODEsStoreNames        % Names for the storesFluxNames         % Names for the fluxesFluxGroups        % Grouping of fluxes (useful for water balance and output)StoreSigns        % Signs to give to stores (-1 is a deficit store), assumes all 1 if not given% attributes set at the beginning of the simulation% directly by the usertheta             % Set of parametersdelta_t           % time stepS0                % initial store valuesinput_climate     % vector of input climatesolver_opts       % options for numerical solving of ODEs% automatically, based on parameter setstore_min         % store minimum valuesstore_max         % store maximum values% attributes created and updated automatically throughout a% simulationt                 % current timestepfluxes            % vector of all fluxesstores            % vector of all storesuhs               % unit hydrographs and still-to-flow fluxessolver_data       % step-by-step info of solver used and residualsstatus            % 0 = model created, 1 = simulation endedendmethods% This will run as soon as any model object is createdfunction [obj] = MARRMoT_model()obj.isOctave = exist('OCTAVE_VERSION', 'builtin')~=0;% if running in Octave, load the optim packageif obj.isOctave; pkg load optim; endendfunction [] = set.isOctave(obj, value); obj.isOctave = value; end% Set methods with checks on inputs for attributes set by the user:function [] = set.delta_t(obj, value)if numel(value) == 1 || isempty(value)obj.delta_t = value;obj.reset();elseerror('delta_t must be a scalar')endendfunction [] = set.theta(obj, value)if numel(value) == obj.numParams || isempty(value)obj.theta = value(:);obj.reset();elseerror(['theta must have ' int2str(obj.numParams) ' elements'])endendfunction [] = set.input_climate(obj, value)if isstruct(value)if isfield(value, 'delta_t')obj.delta_t = value.delta_t;elseif isempty(obj.delta_t)error(['delta_t is not in input climate struct: '...'add it with obj.delta_t = delta_t'])endif isfield(value, {'precip' 'pet' 'temp'})P = value.precip/obj.delta_t;Ea = value.pet/obj.delta_t;T = value.temp/obj.delta_t;obj.input_climate = [P(:) Ea(:) T(:)];obj.reset();elseerror(['Input climate struct must contain fields: '...'precip, pet, temp']);endelseif isnumeric(value)if size(value,2)  || isempty(value)obj.input_climate = value;obj.reset();elseerror(['Input climate must have 3 columns: '...'precip, pet, temp']);endelseerror(['Input climate must either be a struct or '...'a numeric array of 3 columns']);endendfunction [] = set.S0(obj, value)if numel(value) == obj.numStores || isempty(value)obj.S0 = value(:);obj.reset();elseerror(['S0 must have ' int2str(obj.numStores) ' elements'])endendfunction [] = set.solver_opts(obj, value)% add options to default onesobj.solver_opts = obj.add_to_def_opts(value);obj.reset();end% INIT_ runs before each model run to initialise store limits,% auxiliary parameters etc. it calls INIT which is model specificfunction obj = init_(obj)% min and max of storesobj.store_min = zeros(obj.numStores,1);obj.store_max = inf(obj.numStores,1);% empty vectors of fluxes and storest_end = size(obj.input_climate, 1);obj.stores = zeros(t_end, obj.numStores);obj.fluxes = zeros(t_end, obj.numFluxes);% empty struct with the solver dataobj.solver_data.resnorm   = zeros(t_end,1);obj.solver_data.solver = zeros(t_end,1);if(~obj.isOctave); obj.solver_data.solver = categorical(obj.solver_data.solver); end;obj.solver_data.iter   = zeros(t_end,1);% model specific initialisationobj.init();end% RESET is called any time that a user-specified input is changed% (t, delta_t, input_climate, S0, solver_options) and resets any% previous simulation ran on the object.% This is to prevent human error in analysing results.function obj = reset(obj)obj.t = [];                 % current timestepobj.fluxes = [];            % vector of all fluxesobj.stores = [];            % vector of all storesobj.uhs = [];               % unit hydrographs and still-to-flow fluxesobj.solver_data = [];       % step-by-step info of solver used and residualsobj.status = 0;             % 0 = model created, 1 = simulation endedend% ODE approximation with Implicit Euler time-stepping schemefunction err = ODE_approx_IE(obj, S)S = S(:);delta_S = obj.model_fun(S);if obj.t == 1; Sold = obj.S0(:);else; Sold = obj.stores(obj.t-1,:)';enderr = (S - Sold)/obj.delta_t - delta_S';end % SOLVE_STORES solves the stores ODEs function [Snew, resnorm, solver, iter] = solve_stores(obj, Sold)solver_opts = obj.solver_opts;% This reduces the tolerance to a fraction of the smallest store,% if stores are very small, with 1E-6 as minimum% (if resnorm_tolerance is 0.1 as default)resnorm_tolerance = solver_opts.resnorm_tolerance * min(min(abs(Sold)) + 1E-5, 1);% create vectors for each of the three solutions (NewtonRaphon,% fsolve and lsqnonlin), this way if all three run it takes the% best at the end and not the last one.Snew_v    = zeros(3, obj.numStores);resnorm_v = Inf(3, 1);iter_v    = ones(3,1);% first try to solve the ODEs using NewtonRaphsonif(~obj.isOctave) %if MATLAB[tmp_Snew, tmp_fval] = ...NewtonRaphson(@obj.ODE_approx_IE,...Sold,...solver_opts.NewtonRaphson);else              % if Octave[tmp_Snew, tmp_fval] = ...NewtonRaphson_octave(@obj.ODE_approx_IE,...Sold,...solver_opts.NewtonRaphson);endtmp_resnorm = sum(tmp_fval.^2);Snew_v(1,:)  = tmp_Snew;resnorm_v(1) = tmp_resnorm;% if NewtonRaphson doesn't find a good enough solution, run FSOLVEif tmp_resnorm > resnorm_tolerance  [tmp_Snew,tmp_fval,~,tmp_iter] = ...obj.rerunSolver('fsolve', ...              tmp_Snew, ...                  % recent estimatesSold);                         % storages at previous time steptmp_resnorm = sum(tmp_fval.^2);Snew_v(2,:)  = tmp_Snew;resnorm_v(2) = tmp_resnorm;iter_v(2)    = tmp_iter;% if FSOLVE doesn't find a good enough solution, run LSQNONLINif tmp_resnorm > resnorm_tolerance[tmp_Snew,tmp_fval,~,tmp_iter] = ...obj.rerunSolver('lsqnonlin', ...              tmp_Snew, ...                  % recent estimatesSold);                         % storages at previous time steptmp_resnorm = sum(tmp_fval.^2);Snew_v(3,:)  = tmp_Snew;resnorm_v(3) = tmp_resnorm;iter_v(3)    = tmp_iter;endend% get the best solution[resnorm, solver_id] = min(resnorm_v);Snew = Snew_v(solver_id,:);iter = iter_v(solver_id);if(obj.isOctave)solver = solver_id;elsesolvers = ["NewtonRaphson", "fsolve", "lsqnonlin"];solver = solvers(solver_id);endend% RERUNSOLVER Restarts a root-finding solver with different % starting pointsfunction [ Snew, fval, stopflag, stopiter ] = ...rerunSolver( obj,...solverName,...initGuess,...Sold)% get out useful attributessolver_opts = obj.solver_opts.(solverName);solve_fun = @obj.ODE_approx_IE;max_iter = obj.solver_opts.resnorm_maxiter;resnorm_tolerance = obj.solver_opts.resnorm_tolerance * min(min(abs(Sold)) + 1E-5, 1);% Initialize iteration counter, sampling checker and find number of ODEsiter      = 1;resnorm   = resnorm_tolerance + 1;                                 % i.e. greater than the required accuracynumStores = obj.numStores;stopflag  = 1;                                                     % normal function run% Initialise vector of Snew and fval for each iteration, this way you can% keep the best one, not the last one.Snew_v    = zeros(numStores, max_iter);fval_v    = inf(numStores,max_iter);resnorm_v = inf(1, max_iter);Snew = -1 * ones(numStores, 1);% Start the re-sampling% Re-sampling uses different starting points for the solver to see if% solution accuracy improves. Starting points are alternated as follows:% 1. location where the solver got stuck% 2. storages at previous time step% 3. minimum values% 4. maximum values% 5. randomized values close to solution of previous time stepswhile resnorm > resnorm_tolerance% Select the starting pointsswitch itercase 1x0 = initGuess(:);                                     % 1. Location where solver got stuckcase 2x0 = Sold(:);                                          % 2. Stores at t-1case 3x0 = max(-2*10^4.*ones(numStores,1),obj.store_min(:)); % 3. Low values (store minima or -2E4)case 4x0 = min(2*10^4.*ones(numStores,1),obj.store_max(:));  % 4. High values (store maxima or 2E4)otherwisex0 = max(zeros(numStores,1),...Sold(:)+randn(numStores,1).*Sold(:)/10);      % 5. Randomized values close to starting locationend% Re-run the solverif strcmpi(solverName, 'fsolve')[Snew_v(:,iter), fval_v(:,iter), stopflag] = ...fsolve(solve_fun, x0, solver_opts);elseif strcmpi(solverName, 'lsqnonlin')[Snew_v(:,iter), ~,  fval_v(:,iter), stopflag] = ...lsqnonlin(solve_fun, x0, obj.store_min, [], solver_opts);elseerror('Only fsolve and lsqnonlin are supported');endresnorm_v(iter) = sum(fval_v(:,iter).^2);[resnorm,stopiter] = min(resnorm_v);fval = fval_v(:,stopiter);Snew = Snew_v(:,stopiter);% Break out of the loop of iterations exceed the specified maximumif iter >= max_iterstopflag = 0;                                                          % function stopped due to iteration countbreakend% Increase the iteration counteriter = iter + 1;endend% RUN runs the model with a given climate input, initial stores,% parameter set and solver settings.% none of the arguments are needed, they can be set beforehand with% obj.theta = theta; obj.input_climate = input_climate; etc.% then simply obj.run() without arguments.function [] = run(obj,...input_climate,...         S0,...theta,...solver_opts)if nargin > 4 && ~isempty(solver_opts)obj.solver_opts = solver_opts;endif nargin > 3 && ~isempty(theta)obj.theta = theta;endif nargin > 2 && ~isempty(S0)obj.S0 = S0;endif nargin > 1 && ~isempty(input_climate)obj.input_climate = input_climate;end% run INIT_ method, this will calculate all auxiliary parameters% and set up routing vectors and store limitsobj.init_();t_end = size(obj.input_climate, 1);for t = 1:t_endobj.t = t;if t == 1; Sold = obj.S0(:);else; Sold = obj.stores(t-1,:)';end[Snew,resnorm,solver,iter] = obj.solve_stores(Sold);[dS, f] = obj.model_fun(Snew);obj.fluxes(t,:) = f * obj.delta_t;obj.stores(t,:) = Sold + dS' * obj.delta_t;obj.solver_data.resnorm(t) = resnorm;obj.solver_data.solver(t) = solver;obj.solver_data.iter(t) = iter;obj.step();endobj.status = 1;end% GET_OUTPUT runs the model exactly like RUN, but output is% consistent with current MARRMoTfunction [fluxOutput,...fluxInternal,...storeInternal,...waterBalance,...solverSteps] = get_output(obj,...varargin)if nargin > 1 || isempty(obj.status) || obj.status == 0 obj.run(varargin{:});end% --- Fluxes leaving the model ---          fg = fieldnames(obj.FluxGroups);fluxOutput = struct();for k=1:numel(fg)idx = abs(obj.FluxGroups.(fg{k}));signs = sign(obj.FluxGroups.(fg{k}));fluxOutput.(fg{k}) = sum(signs.*obj.fluxes(:,idx),2);end% --- Fluxes internal to the model ---fluxInternal = struct;for i = 1:obj.numFluxesfluxInternal.(char(obj.FluxNames{i})) = obj.fluxes(:,i)';end% --- Stores ---storeInternal = struct;for i = 1:obj.numStoresstoreInternal.(char(obj.StoreNames{i})) = obj.stores(:,i)';end% --- Water balance, if requested ---if nargout >= 4waterBalance = obj.check_waterbalance();end% --- step-by-step data of the solver, if requested ---if nargout == 5solverSteps = obj.solver_data;endend% CHECK_WATERBALANCE returns the waterbalance% like in MARRMoT1, it will print to screenfunction [out] = check_waterbalance(obj, varargin)if nargin > 1 || isempty(obj.status) || obj.status == 0 obj.run(varargin{:});end% Get variablesP  = obj.input_climate(:,1);fg = fieldnames(obj.FluxGroups);OutFluxes = zeros(1,numel(fg));for k=1:numel(fg)                                              % cumulative of each flow leaving the modelidx = abs(obj.FluxGroups.(fg{k}));signs = sign(obj.FluxGroups.(fg{k}));OutFluxes(k) = sum(sum(signs.*obj.fluxes(:,idx), 1),2);endif isempty(obj.StoreSigns); obj.StoreSigns = repelem(1, obj.numStores); enddS = obj.StoreSigns(:) .* (obj.stores(end,:)' - obj.S0);          % difference of final and initial storage for each storeif isempty(obj.uhs); obj.uhs = {}; endR = cellfun(@(uh) sum(uh(2,:)), obj.uhs);                      % cumulative of each flows still to be routed% calculate water balanceout = sum(P) - ...                                             % input from precipitationsum(OutFluxes) - ...                                       % all fluxes leaving the model (some may be entering, but they should have a negative sign)sum(dS) - ...                                              % all differences in storagesum(R);                                                    % all flows still being routeddisp(['Total P  = ',num2str(sum(P)),' mm.'])for k = 1:numel(fg)disp(['Total ',char(fg(k)),' = ',...num2str(-OutFluxes(k)),' mm.'])endfor s = 1:obj.numStoresif obj.StoreSigns(s) == -1ending=' (deficit store).';elseending='.';enddisp(['Delta S',num2str(s),' = ',...num2str(-dS(s)),' mm',ending])endif ~isempty(R)disp(['On route = ',num2str(-sum(R)),' mm.'])enddisp('-------------')disp(['Water balance = ', num2str(out), ' mm.'])end% GET_STREAMFLOW only returns the streamflow, runs the model if it% hadn't run already.function Q = get_streamflow(obj, varargin)if nargin > 1 || isempty(obj.status) || obj.status == 0obj.run(varargin{:});endQ = sum(obj.fluxes(:,obj.FluxGroups.Q),2);end% CALIBRATE uses the chosen algorithm to find the optimal parameter% set, given model inputs, objective function and observed streamflow.% the function chosen in algorithm should have the same inputs and% outputs as MATLAB's fminsearch.function  [par_opt,...                                             % optimal parameter setof_cal,...                                              % value of objective function at par_optstopflag,...                                            % flag indicating reason the algorithm stoppedoutput] = ...                                           % output, see fminsearch for detailcalibrate(obj,...Q_obs,...                           % observed streamflowcal_idx,...                         % timesteps to use for model calibrationoptim_fun,...                       % function to use for optimisation (must have same structure as fminsearch)par_ini,...                         % initial parameter estimatesoptim_opts,...                      % options to optim_funof_name,...                         % name of objective function to useinverse_flag,...                    % should the OF be inversed?display,...                         % should I display information about the calibration?varargin)                           % additional arguments to the objective functionif isempty(obj.input_climate) || isempty(obj.delta_t) ||...isempty(obj.S0) || isempty(obj.solver_opts)error(['input_climate, delta_t, S0 and solver_opts '...'attributes must be specified before calling '...'calibrate.']);end% if the list of timesteps to use for calibration is empty,% use all steps if isempty(cal_idx)cal_idx = 1:length(Q_obs);end% use display by defaultif isempty(display)display = true;end% use the data from the start to the last value of cal_idx to% run the simulationif islogical(cal_idx); cal_idx = find(cal_idx); endinput_climate_all = obj.input_climate;obj.input_climate = input_climate_all(1:max(cal_idx),:);Q_obs = Q_obs(1:max(cal_idx));% if the initial parameter set isn't set,  start from mean% values of parameter rangeif isempty(par_ini)par_ini = mean(obj.parRanges,2);end% helper function to calculate fitness given a set of% parametersfunction fitness = fitness_fun(par)Q_sim = obj.get_streamflow([],[],par);fitness = (-1)^inverse_flag*feval(of_name, Q_obs, Q_sim, cal_idx, varargin{:});end% display some useful things for the user to make sure they% used the right settingsif displaydisp('---')disp(['Starting calibration of model ' class(obj) '.'])disp(['Simulation will run for timesteps 1-' num2str(max(cal_idx)) '.'])% this is a bit ugly, but it formats the list of cal_idx in% a pretty and concise waycal_idx = sort(cal_idx);i = 1;previous = cal_idx(i);cal_idx_str = num2str(previous);  while i < numel(cal_idx)  i = i + 1;  if cal_idx(i)-previous == 1  i = find(diff(cal_idx(i:end)) ~= 1, 1) + i - 1;  if isempty(i); i = numel(cal_idx); end  previous = cal_idx(i);  cal_idx_str = [cal_idx_str '-', num2str(previous)];  else  previous = cal_idx(i);  cal_idx_str = [cal_idx_str ', ', num2str(previous)];  end  enddisp(['Objective function ' of_name ' will be calculated in time steps ' cal_idx_str '.'])disp(['The optimiser ' optim_fun ' will be used to optimise the objective function.'])disp(['Options passed to the optimiser:'])disp(optim_opts)disp('All other options are left to their default,')disp('check the source code of the optimiser to find these default values.')disp('---')end[par_opt,...                                                  % optimal parameter set at the end of the optimisationof_cal,...                                                % value of the objective function at par_optstopflag,...                                              % flag indicating reason the algorithm stoppedoutput] = ...                                             % output, see fminsearch for detailfeval(optim_fun,...                             % run the optimisation algorithm chosen@fitness_fun,...                          % function to optimise is the fitness functionpar_ini,...                               % initial parameter setoptim_opts);                              % optimiser options% if of_cal was inverted, invert it back before returningof_cal = (-1)^inverse_flag * of_cal;% reset the whole input climate as it was before the% calibrationobj.input_climate = input_climate_all;end% function to return default solver optionsfunction solver_opts = default_solver_opts(obj)solver_opts.resnorm_tolerance = 0.1;                                       % Root-finding convergence tolerancesolver_opts.resnorm_maxiter   = 6;                                         % Maximum number of re-runs used in rerunSolversolver_opts.NewtonRaphson = optimset('MaxIter', obj.numStores * 10);% if MATLABif(~obj.isOctave)solver_opts.fsolve = optimoptions('fsolve',...'Display','none',...                     % Disable display settings'JacobPattern', obj.JacobPattern);solver_opts.lsqnonlin = optimoptions('lsqnonlin',...                       % lsqnonlin settings for cases where fsolve fails'Display','none',...'JacobPattern',obj.JacobPattern,...'MaxFunEvals',1000);else % if OCTAVEsolver_opts.fsolve = optimset('Display', 'off');solver_opts.lsqnonlin = optimset('Display', 'off', 'MaxFunEvals',1000);endend% function to add new solver opts to the default onesfunction solver_opts = add_to_def_opts(obj, opts)def_opts = obj.default_solver_opts();             if nargin == 1 || isempty(opts)solver_opts = def_opts;elsedef_fields = fieldnames(def_opts);% for each field in the default options (5 at the moment)for k = 1:length(def_fields)field = def_fields{k};% if the field is not provided, use the default oneif ~isfield(opts, field) || isempty(opts.(field))solver_opts.(field) = def_opts.(field);% if the field is provided, and the dafault is a struct,% add the new values to the default structelseif isstruct(def_opts.(field))solver_opts.(field) = optimset(def_opts.(field),opts.(field));% if the field is provided, and the dafault is not a struct,% discard the default and use the new valueelsesolver_opts.(field) = opts.(field);endendend             endend
end