python编程练习：模拟水文模型中的水箱模型(tank model)，不含参数率定过程

一、水箱模型结构

二、代码

import matplotlib.pyplot as plt
import numpy as np
import math
import time
# 设置标签为负号可显示，坐标轴可显示中文
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False# 实测数据，time为时刻，rain为降雨强度，qoc为产流
ori_time = ['0:00', '1:00', '2:00', '3:00', '4:00', '5:00', '6:00', '7:00', '8:00', '9:00', '10:00', '11:00', '12:00']
ori_rain = [5, 5, 5, 5, 5, 5, 5, 5, 5, 3, 2, 1, 0]
ori_qoc = [0.477, 0.563, 0.744, 1.18, 1.966, 3.315, 4.269, 5.271, 5.859, 4.943, 4.365, 3.97, 3.758]def cre_single_tank(in_rain, d_factor, r_factor, thres, h_intial=None):""":param h_intial: 单孔水箱的初始贮水量，默认为10:param h_initial: 单个水箱的初始存储量,默认为10:param in_rain: 实测的rainfall数据:param d_factor: 单侧孔水箱的下渗系数:param r_factor: 单侧孔水箱的流出系数:param thres: 单侧孔水箱的孔高:return: out_down为下渗量的时间变化；out_right为流出量的时间变化"""if h_intial is None:h_intial = 0out_down = []out_right = []z_end = []# 计算初始时的z_endtemp = in_rain[0] + h_intialout_down.append(temp * d_factor)if temp <= thres:out_right.append(0)else:out_right.append((temp - thres) * r_factor)z_end.append(temp - out_down[0] - out_right[0])# 计算其他时刻的z_endfor i in range(1, len(in_rain)):z_start = in_rain[i] + z_end[i-1]down = z_start * d_factorif z_start <= thres:right = 0else:right = (z_start - thres) * r_factorz_end.append(z_start - down - right)out_down.append(down)out_right.append(right)out_down = np.array([round(i, 3) for i in out_down])out_right = np.array([round(i, 3) for i in out_right])return out_down, out_rightdef cre_double_tank(in_rain,  d_factor, r1_factor, r2_factor, thres_1, thres_2, h_intial=None):""":param h_initial: 双孔水箱的初始存储量,默认为10:param in_rain: 实测的rainfall数据:param d_factor: 双侧孔水箱的下渗系数:param r1_factor: 双侧孔水箱下孔的流出系数:param r2_factor: 双侧孔水箱上孔的流出系数:param thres_1: 双侧孔水箱下孔的孔高:param thres_2: 双侧孔水箱上孔的孔高:return: out_down为下渗量的时间变化；out_right为流出量的时间变化"""out_down = []out_right = []z_end = []if h_intial is None:h_intial = 0# 计算初始时的z_endtemp = in_rain[0] + h_intialif temp <= thres_1:right_1 = 0else:right_1 = (temp - thres_1) * r1_factorif temp <= thres_2:right_2 = 0else:right_2 = (temp - thres_2) * r2_factorout_down.append(temp * d_factor)out_right.append(right_1 + right_2)z_end.append(temp - out_down[0] - out_right[0])# 计算其他时刻的z_endfor i in range(1, len(in_rain)):z_start = in_rain[i] + z_end[i - 1]down = z_start * d_factorif z_start <= thres_1:right_1 = 0else:right_1 = (z_start - thres_1) * r1_factorif z_start <= thres_2:right_2 = 0else:right_2 = (z_start - thres_2) * r2_factorz_end.append(z_start - down - right_1 - right_2)out_down.append(down)out_right.append(right_1 + right_2)out_down = np.array([round(i, 3) for i in out_down])out_right = np.array([round(i, 3) for i in out_right])return out_down, out_rightdef cal_correlation(real_data, sim_data):x_bar = np.mean(real_data)y_bar = np.mean(sim_data)ssr = 0var_x = 0var_y = 0for i in range(len(real_data)):diff_xx_bar = real_data[i] - x_bardiff_yy_bar = sim_data[i] - y_barssr += (diff_xx_bar * diff_yy_bar)var_x += diff_xx_bar ** 2var_y += diff_yy_bar ** 2sst = math.sqrt(var_x * var_y)return round(ssr/sst, 4)def cal_dy(real_data, sim_data):num_data = len(real_data)aver_real = sum(real_data) / num_datat_1 = 0t_2 = 0for i in range(num_data):t_1 += (real_data[i] - sim_data[i]) ** 2t_2 += (real_data[i] - aver_real) ** 2r = 1 - (t_1 / t_2)return round(r, 2)def do_simulation(switchs=None):# F_n: 第n层水箱下渗量的时间序列# R_n: 第n层水箱流出量的时间序列if switchs is None:switchs = [1, 0]# 搭建水箱模型F_1, R_1 = cre_double_tank(ori_rain, d_factors[0], r_factors[0], r_factors[1], hole_heights[0], hole_heights[1], h_start[0])F_2, R_2 = cre_single_tank(F_1, d_factors[1], r_factors[2], hole_heights[2], h_start[1])F_3, R_3 = cre_single_tank(F_2, d_factors[2], r_factors[3], hole_heights[3], h_start[2])F_4, R_4 = cre_single_tank(F_3, d_factors[3], r_factors[4], hole_heights[4], h_start[3])# cal_qoc: 模拟结果得到qoc时间序列cal_qoc = R_1 + R_2 + R_3 + R_4cal_qoc = [round(qoc, 3) for qoc in cal_qoc]# 用R-square和确定性系数dy衡量模拟效果的好坏r_square.append(cal_correlation(ori_qoc, cal_qoc))dy.append(cal_dy(ori_qoc, cal_qoc))# 进行绘图figure_show_switch = switchs[0]types = ['kx-', 'ks-']each_layer_switch = 1if figure_show_switch:plt.plot(ori_qoc, types[0], linewidth=1, label='实测', markeredgecolor='blue')plt.plot(cal_qoc, types[1], linewidth=1, label='模拟')if each_layer_switch:plt.plot(R_1, 'r', linewidth=1, label='L1')plt.plot(R_2, 'g', linewidth=1, label='L2')plt.plot(R_3, 'b', linewidth=1, label='L3')plt.plot(R_4, 'y', linewidth=1, label='L4')plt.legend()plt.show()# 生成模拟日志save_log_switch = switchs[1]time_run = time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(time.time()))log_name = 'log.txt'if save_log_switch:with open(log_name, 'a+') as f_obj:f_obj.writelines('{0}\n'.format("=" * 30))# 当前时间f_obj.writelines(time_run + '\n')f_obj.writelines("\nInput parameters:\n")f_obj.writelines("d_factors = {0}\n".format(d_factors))f_obj.writelines("r_factors = {0}\n".format(r_factors))f_obj.writelines("hole_heights = {0}\n".format(hole_heights))f_obj.writelines("\nSimulated result:\n")f_obj.writelines("cal_qoc = {0}\n".format(list(cal_qoc)))f_obj.writelines("R_square:{0}\n".format(r_square[-1]))f_obj.writelines("dy:{0}\n".format(dy[-1]))f_obj.writelines("=" * 30)# 输出拟合结果print("r_square: {0}".format(r_square[-1]))print("dy: {0}".format(dy[-1]))return max(dy)"""
与四层水箱模型相关的所有参数如下孔高         下渗系数     流出系数     初始贮水量
第一层  hole_h[0]/hole_h[1]   d_f[0]   r_f[0]/r_f[1]  h_start[0]
第二层        hole_h[2]       d_f[1]       r_f[2]     h_start[1]
第三层        hole_h[3]       d_f[2]       r_f[3]     h_start[2]
第四层        hole_h[4]       d_f[3]       r_f[4]     h_start[3]
"""
# 默认的取值
d_factors = [0.2, 0.1, 0.05, 0]
r_factors = [0.2, 0.1, 0.05, 0.03, 0.01]
hole_heights = [40, 25, 10, 15, 15]
h_start = [0, 0, 0, 0]
r_square = []
dy = []# 第一层参数
hole_heights[0] = 45
hole_heights[1] = 15
r_factors[0] = 0.3
r_factors[1] = 0.06
d_factors[0] = 0.2# 第二层参数
hole_heights[2] = 10
r_factors[2] = 0.06
d_factors[1] = 0.1# 第三层参数
hole_heights[3] = 0
r_factors[3] = 0.03
d_factors[2] = 0.03# 第四层参数
hole_heights[4] = 0
r_factors[4] = 0.005
d_factors[3] = 0# 进行模拟，两个1分别控制绘图与日志
do_simulation([1, 1])

三、模拟结果
可以输出模拟后的产流序列与实测产流序列的对比图，L1~L4为模拟时各层水箱对产流序列。

生成日志文件如下图

接着通过手动率定水箱模型中的各个参数，使确定系数dy(-1~1)尽可能接近1，使拟合程度增加至符合要求。

python编程练习：模拟水文模型中的水箱模型(tank model)，不含参数率定过程相关推荐

NCT青少年编程能力等级测试Python编程二级-模拟卷（含答案）
参考答案在文章后边部分,请看到后半部分的答案分割线,非常感谢哦! 试题NCT-Python编程三级-模拟卷1(含答案一.选择题 1．Python语言属于( ). A．机器语言 B．汇编语言 C． ...
NCT青少年编程能力等级测试Python编程三级-模拟卷1（含答案）
参考答案在文章后边部分,请看到后半部分的答案分割线,非常感谢哦! 试题NCT-Python编程三级-模拟卷2(含答案练习一.选择题 1．下面( )是Python合法的变量名 A．int32 B ...
NCT青少年编程能力等级测试Python编程一级-模拟卷（含答案）
参考答案在文章后边部分,请看到后半部分的答案分割线,非常感谢哦! 一.程序填空 1．编写程序,实现从键盘输入数据,数据前三位的ASCII值加2,从第四位开始ASCII值加3. 2．请在空格处填写正确的 ...
python编程字典100例_python中字典(Dictionary)用法实例详解
本文实例讲述了python中字典(Dictionary)用法.分享给大家供大家参考.具体分析如下: 字典(Dictionary)是一种映射结构的数据类型,由无序的"键-值对"组成. ...
如何在matlab中建立水箱模型_在MATLAB中实现水箱液位控制系统的设计
在 MATLAB 中实现水箱液位控制系统的设计 [摘要] 本论文的目的是设计双容水箱液位串级控制系统. 在设计中充分利用计算机技术, 自动控制技术, 以实现对水箱液位的串级控制. 首先对被控对象的 ...
python编程手机模拟点击_python模拟点击玩游戏的实例讲解
小编发现很多小伙伴都喜欢玩一些游戏,而手游因为玩的场景限制不多,所以受众的人更多.游戏里有很多重复的任务需要我们完成,虽然过程非常无聊,但是为了任务奖励还是有很多小伙伴不厌其烦的去做.那么,有没有什么 ...
ipad python编程软件_在iPad中运行Python
从一个喜欢编程的人的角度看,任何移动设备其实都只是"可编程计算器". iPad上也有一个Python的解释器的应用Python for iOS.不过,让人气愤的是,这个应用竟然是收 ...
python编程手机模拟点击_python简单的模拟点击（一）
使用python模拟百度搜索"data_bug"的博客下面根据代码一步步带你解释 from selenium import webdriver from selenium.web ...
linux环境下python编程指南,在Linux系统中搭建Python编程环境
Linux系统是为编程而设计的,因此在大多数Linux计算机中都默认安装了Python. 1. 检查Python版本在系统中运行应用程序Terminal(如果是Ubuntu,可按Ctrl+Alt+T ...

python编程练习：模拟水文模型中的水箱模型(tank model)，不含参数率定过程

python编程练习：模拟水文模型中的水箱模型(tank model)，不含参数率定过程相关推荐

最新文章

热门文章