import pandas as pd
import numpy as np
import os
import time
#1.参数设置与数据输入
data = pd.read_excel(r"C:\Users\80401\Desktop\SHUJU2.xlsx")
P = data['P'].values # 转换为 NumPy 数组
E0 = data['E0'].values
# 参数初始化
K = float(1.20)
C = float(0.16)
Im = float(0.003)
WM = float(140)
WUM = float(20)
WLM = float(60)
WDM = float(60)
b = float(0.3)
WU = float(10.0)
WL = float(40.0)
WD = float(60.0)
# 预处理:划分不透水面积和透水面积
P_im = P * Im # 不透水面积降雨
P_soil = P * (1 - Im) # 透水面积降雨
# 结果存储列表
results = []
#2.三层蒸发计算模块
def evp(WU, WL, WD, Ep, P):
"""三层蒸发计算模块"""
# 初始化所有变量
EU = EL = ED = 0.0
if P == 0 and WU == 0:
# 上层无水时,考虑下层蒸发
if WL >= C * WLM:
EL = min(Ep, WL) # 限制蒸发量不超过可用水量
elif WL >= C * Ep:
EL = C * Ep
else:
EL = WL
remaining_ep = max(Ep - EL, 0) # 确保蒸发需求非负
ED = min(C * remaining_ep, WD)
elif WU + P >= Ep: # 三层蒸散发
EU = Ep
else:
EU = WU + P
remaining_ep = Ep - EU
if WL >= C * WLM:
EL = min(remaining_ep * WL / WLM, WL)
elif WL >= C * remaining_ep:
EL = C * remaining_ep
else:
EL = WL
ED = min(C * (remaining_ep - EL), WD)
E = EU + EL + ED
PE = P - E
return EU, EL, ED, E, PE
#3.蓄满产流计算模块
def runoff(PE,W,WM,b):
"""蓄满产流计算"""
if PE <= 0:
return 0.0
# 完全饱和时所有净雨转化为径流
if W >= WM:
return PE
WMM = WM * (1 + b)
a = WMM * (1 - (1 -W / WM) ** (1 / (1 + b)))
if a + PE <= WMM:
R = PE + W - WM + WM * (1 - (PE + a) / WMM) ** (1 + b)
else:
R = PE + W - WM
return max(R, 0) # 确保产流不为负
#4.土壤蓄水量更新模块
def update(WU, WL, WD, EU, EL, ED, PE, R):
"""更新土壤含水量"""
# 1. 计算实际入渗量(净雨量扣除产流)
net_water = PE - R
# 2. 分层补充水量 (先补充上层,再下层,最后深层)
# 上层
add_WU = min(WUM - WU, max(0, net_water))
new_WU = max(0, min(WU + add_WU, WUM))
net_water -= add_WU #下渗量
# 下层
add_WL = min(WLM - WL, max(0, net_water))
new_WL = max(0, min(WL + add_WL, WLM))
net_water -= add_WL #下渗量
# 深层
add_WD = min(WDM - WD, max(0, net_water))
new_WD = max(0, min(WD + add_WD, WDM))
# 3. 水量平衡检查 (使用透水面积部分)
total_change = (new_WU - WU) + (new_WL - WL) + (new_WD - WD)
expected_change = PE - R # 净雨量-产流
# 5.调整误差(确保水量平衡)
if abs(total_change - expected_change) > 0.001:
adjustment = expected_change - total_change
# 将调整量加到深层(或按需分配到各层)
new_WD = max(0, min(new_WD + adjustment, WDM))
return new_WU, new_WL, new_WD
# 主循环
for i in range(len(P)):
try:
# 获取当前时段的原始降雨和蒸发能力
current_P_total = P[i] # 原始总降雨量
current_E0 = E0[i]
# 划分不透水面积和透水面积
current_P_im = P_im[i] # 不透水面积降雨
current_P_soil = P_soil[i] # 透水面积降雨
# 计算实际蒸发能力
current_Ep = current_E0 * K
# 三层蒸散发计算 (使用透水面积降雨)
EU, EL, ED, E, PE = evp(WU, WL, WD, current_Ep, current_P_soil)
# 产流计算 (仅透水面积)
W = WU + WL + WD
R_soil = runoff(PE, W, WM, b)
# 总产流 = 透水面积产流 + 不透水面积直接径流
R_total = R_soil + current_P_im
# 更新土壤含水量 (仅透水面积部分)
new_WU, new_WL, new_WD = update(WU, WL, WD, EU, EL, ED, PE, R_soil)
# 记录结果
results.append({
'Period': i + 1,
'P_total': round(current_P_total, 3),
'P_im': round(current_P_im, 3),
'P_soil': round(current_P_soil, 3),
'E0': round(current_E0, 2),
'Ep': round(current_Ep, 2),
'EU': round(EU, 3),
'EL': round(EL, 3),
'ED': round(ED, 3),
'E': round(E, 3),
'PE': round(PE, 3),
'R_soil': round(R_soil, 3),
'R_im': round(current_P_im, 3),
'R_total': round(R_total, 3),
'WU': round(new_WU, 3),
'WL': round(new_WL, 3),
'WD': round(new_WD, 3),
'Total_W': round(new_WU + new_WL + new_WD, 3)
})
# 更新状态变量
WU, WL, WD = new_WU, new_WL, new_WD
except Exception as e:
print(f"错误发生在第 {i + 1} 时段: {str(e)}")
print(f"当前参数: WU={WU}, WL={WL}, WD={WD}, Ep={Ep[i]}, P={P[i]}")
break
# 转换为DataFrame并保存结果
results_df = pd.DataFrame(results)
print("\n计算结果:")
print(results_df)
# 保存到Excel
results_df.to_excel("水文模型计算结果KESHE.xlsx", index=False)
print("结果已保存到 '水文模型计算结果KESHE.xlsx'")
# # 打印前10个时段的结果
# print("\n前10个时段计算结果:")
# print(results_df.head(10))
#
# # 保存到桌面(避免权限问题)
# desktop_path = os.path.join(os.path.expanduser("~"), "Desktop")
# output_file = os.path.join(desktop_path, f"水文模型计算结果_{time.strftime('%Y%m%d%H%M')}.xlsx")
#
# # 使用openpyxl引擎确保写入成功
# results_df.to_excel(output_file, index=False, engine='openpyxl')
# print(f"\n结果已保存到: {output_file}")
#
# # 打印汇总统计
# print("\n模型运行汇总:")
# print(f"总时段数: {len(results)}")
# print(f"总降雨量: {round(sum(P), 2)} mm")
# print(f"总蒸发量: {round(sum(results_df['E']), 2)} mm")
# print(f"总产流量: {round(sum(results_df['R']), 2)} mm")
最新发布