import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# ====================== 参数设置 ======================
# 流域参数
AREA = 553 # 流域面积(km²)
DT = 3 # 计算时段长度(小时)
# 产流参数 (示例值，实际应通过任务1优选)
params = {
'WM': 140, # 张力水蓄水容量(mm)
'WUM': 20, # 上层张力水容量(mm)
'WLM': 60, # 下层张力水容量(mm)
'WDM': 60, # 深层张力水容量(mm)
'b': 0.3, # 张力水蓄水容量曲线指数
'c': 0.16, # 深层蒸发系数
'Kc': 0.85, # 蒸发折算系数(应通过优化确定)
'fc': 10.0, # 稳定下渗率(mm/3h)
'WU0': 20, # 初始上层张力水(mm) - 取容量值
'WL0': 60, # 初始下层张力水(mm) - 取容量值
'WD0': 60 # 初始深层张力水(mm) - 取容量值
}
# 汇流参数
unit_hydrograph = [0, 40, 80, 130, 100, 80, 48, 20, 10, 5, 0] # 单位线(m³/s)
CG = 0.968 # 地下径流出流系数
Qg0 = 55.3 # 初始地下径流(m³/s)
# ====================== 新安江模型核心类 ======================
class XinanjiangModel:
def __init__(self, params):
"""初始化模型参数"""
self.WM = params['WM']
self.WUM = params['WUM']
self.WLM = params['WLM']
self.WDM = params['WDM']
self.b = params['b']
self.c = params['c']
self.Kc = params['Kc']
self.fc = params['fc']
# 初始化土壤含水量
self.WU = params['WU0']
self.WL = params['WL0']
self.WD = params['WD0']
def evapotranspiration(self, P, E0):
"""三层蒸发计算"""
Ep = E0 * self.Kc # 蒸散发能力
# 上层蒸发
EU = min(self.WU + P, Ep)
# 剩余蒸发能力
E_remain = Ep - EU
# 下层蒸发
if self.WL > self.c * self.WLM:
EL = min(E_remain, self.WL)
else:
EL = min(self.c * E_remain, self.WL)
# 深层蒸发
E_remain -= EL
ED = min(self.c * E_remain, self.WD)
# 总蒸发量
E_total = EU + EL + ED
# 更新土壤含水量
self.WU = max(0, self.WU + P - EU)
self.WL = max(0, self.WL - EL)
self.WD = max(0, self.WD - ED)
return E_total
def runoff_generation(self, P, E0):
"""产流计算"""
# 计算蒸发
E = self.evapotranspiration(P, E0)
PE = P - E
if PE <= 0:
return 0.0, 0.0, 0.0 # 无产流
# 当前土壤总含水量
W = self.WU + self.WL + self.WD
# 张力水蓄水容量曲线计算
WMM = self.WM * (1 + self.b)
A = WMM * (1 - (1 - min(W / self.WM, 1.0)) ** (1 / (1 + self.b)))
# 产流量计算
if PE + A < WMM:
R = PE - (self.WM - W) + self.WM * (1 - (A + PE) / WMM) ** (1 + self.b)
else:
R = PE - (self.WM - W)
R = max(R, 0) # 确保产流不为负
# 更新土壤含水量
W_new = min(W + PE - R, self.WM)
self.WU = min(W_new * self.WUM / self.WM, self.WUM)
self.WL = min(W_new * self.WLM / self.WM, self.WLM)
self.WD = min(W_new - self.WU - self.WL, self.WDM)
# 分水源计算
if PE > self.fc:
RG = min(R, self.fc) # 地下径流
RD = R - RG # 直接径流
else:
RG = R
RD = 0.0
return R, RD, RG
# ====================== 数据准备 ======================
# 创建DataFrame
data = {
'T': pd.to_datetime([
'2004-9-23 12:00', '2004-9-23 15:00', '2004-9-23 18:00', '2004-9-23 21:00',
'2004-9-24 00:00', '2004-9-24 03:00', '2004-9-24 06:00', '2004-9-24 09:00',
'2004-9-24 12:00', '2004-9-24 15:00', '2004-9-24 18:00', '2004-9-24 21:00',
'2004-9-25 00:00', '2004-9-25 03:00', '2004-9-25 06:00', '2004-9-25 09:00',
'2004-9-25 12:00', '2004-9-25 15:00', '2004-9-25 18:00', '2004-9-25 21:00',
'2004-9-26 00:00', '2004-9-26 03:00', '2004-9-26 06:00', '2004-9-26 09:00',
'2004-9-26 12:00', '2004-9-26 15:00', '2004-9-26 18:00', '2004-9-26 21:00'
]),
'E': [1.3, 1.3, 1.3, 1.3, 1.2, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9,
0.8, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.2, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2],
'P1': [6.2, 7.6, 6.2, 8.8, 25, 29.9, 38.6, 6.9, 28.3, 25.6, 93.9, 85.3,
51.5, 39.8, 43.2, 20.5, 10.5, 7.4, 1.8, 0.2, 0, 0, 0, 0, 0, 0, 0, 0],
'P2': [9.9, 16, 6.4, 17.2, 34.8, 29.2, 24.8, 7.5, 29.9, 42.7, 137.6, 90.8,
47.7, 70.3, 47.3, 13.3, 8, 8.4, 2.8, 0, 0, 0, 0, 0, 0, 0, 0, 0],
'P3': [21.6, 20.6, 14.9, 29.4, 35.3, 43.9, 46.9, 6.1, 34.2, 39.8, 124, 85,
49.2, 42.1, 61.5, 15.8, 1.8, 7.6, 2.1, 0.3, 0, 0, 0, 0, 0, 0, 0, 0],
'P4': [17.3, 12.6, 15.9, 18.5, 24.6, 37.8, 33, 12.3, 28.5, 75.4, 13.2, 75.9,
38.5, 97.7, 45.9, 13.1, 3.3, 10.9, 4.6, 0, 0, 0, 0, 0, 0, 0, 0, 0]
}
df = pd.DataFrame(data)
# 权重系数
weights = [0.33, 0.14, 0.33, 0.20]
# 计算面平均雨量
df['P_avg'] = df[['P1', 'P2', 'P3', 'P4']].dot(weights)
# ====================== 产流计算 ======================
# 初始化模型
model = XinanjiangModel(params)
# 存储结果
df['R'] = 0.0 # 总径流(mm)
df['RD'] = 0.0 # 直接径流(mm)
df['RG'] = 0.0 # 地下径流(mm)
# 逐时段计算产流
for i in range(len(df)):
P = df.loc[i, 'P_avg']
E0_val = df.loc[i, 'E']
R, RD, RG = model.runoff_generation(P, E0_val)
df.loc[i, 'R'] = R
df.loc[i, 'RD'] = RD
df.loc[i, 'RG'] = RG
# ====================== 汇流计算 ======================
# 单位转换系数 (mm/3h → m³/s)
# 公式: 流量(m³/s) = 径流深(mm) * 流域面积(km²) / (3.6 * 时段小时数)
conversion_factor = AREA / (3.6 * DT)
# 1. 直接径流汇流 (单位线法)
# 将径流深转换为单位线倍比
RD_converted = df['RD'].values * conversion_factor / 10 # 单位线为10mm单位线
# 卷积计算
direct_flow = np.convolve(RD_converted, unit_hydrograph, mode='full')[:len(df)]
# 2. 地下径流汇流 (线性水库法)
# 将径流深转换为流量
RG_converted = df['RG'].values * conversion_factor
# 线性水库演算
groundwater_flow = np.zeros(len(df))
groundwater_flow[0] = Qg0 # 初始值
for i in range(1, len(df)):
groundwater_flow[i] = CG * groundwater_flow[i - 1] + (1 - CG) * RG_converted[i]
# 3. 总流量过程
total_flow = direct_flow + groundwater_flow
# ====================== 结果可视化 ======================
plt.figure(figsize=(14, 10))
# 1. 降雨过程
plt.subplot(3, 1, 1)
plt.bar(df['T'], df['P_avg'], width=0.02, color='blue')
plt.title('流域平均降雨过程 (3小时雨量)')
plt.ylabel('降雨量 (mm)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xticks(rotation=45)
# 2. 径流组分
plt.subplot(3, 1, 2)
plt.plot(df['T'], df['R'], 'k-', label='总径流(R)')
plt.plot(df['T'], df['RD'], 'r--', label='直接径流(RD)')
plt.plot(df['T'], df['RG'], 'b-.', label='地下径流(RG)')
plt.title('径流组分过程')
plt.ylabel('径流量 (mm)')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.xticks(rotation=45)
# 3. 流量过程
plt.subplot(3, 1, 3)
plt.plot(df['T'], total_flow, 'k-', linewidth=2, label='总流量')
plt.plot(df['T'], direct_flow, 'r--', label='直接径流')
plt.plot(df['T'], groundwater_flow, 'b-.', label='地下径流')
plt.title('流量过程线')
plt.xlabel('时间')
plt.ylabel('流量 (m³/s)')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.xticks(rotation=45)
plt.tight_layout()
plt.savefig('hydrograph_analysis.png', dpi=300)
plt.show()
# ====================== 结果保存 ======================
# 添加流量结果到DataFrame
df['Direct_Flow'] = direct_flow
df['Groundwater_Flow'] = groundwater_flow
df['Total_Flow'] = total_flow
# 保存结果
df.to_excel(r'C:\Users\80401\Desktop\hydrological_simulation_results.xlsx', index=False)
print("模拟结果已保存到 'hydrological_simulation_results.xlsx'")改进代码，要求不使用scipy库
最新发布
本文探讨了在十亿个数中找到最大的一百万个数的有效算法。介绍了三种方法：直接排序、使用大顶堆及快速选择算法。重点讨论了每种方法的时间复杂度及其适用场景。
扫一扫
885
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
