import pandas as pd
import numpy as np
import os
import matplotlib.pyplot as plt
import time
# ====================== 模型参数设置 ======================
class ModelParams:
"""水文模型参数容器类"""
def __init__(self):
self.Kc = 1.20 # 蒸发系数
self.c = 0.16 # 深层蒸发系数
self.Im = 0.003 # 不透水面积比例
self.WM = 140.0 # 流域蓄水容量
self.WUM = 20.0 # 上层蓄水容量
self.WLM = 60.0 # 中层蓄水容量
self.WDM = 60.0 # 深层蓄水容量
self.b = 0.3 # 蓄水容量曲线指数
self.init_WU = 10.0 # 初始上层含水量
self.init_WL = 40.0 # 初始中层含水量
self.init_WD = 60.0 # 初始深层含水量
# ====================== 核心功能函数 ======================
def calculate_evapotranspiration(P, E0, WU, WL, WD, params):
"""
三层蒸发模型计算
返回: (总蒸发, 上层蒸发, 中层蒸发, 深层蒸发, 上层新含水量, 中层新含水量, 深层新含水量)
"""
# 计算实际蒸散发能力
Ep = E0 * params.Kc
# 1. 上层蒸发计算
if WU + P >= Ep:
EU = Ep
new_WU = WU + P - EU
else:
EU = WU + P
new_WU = 0
remaining_ep = Ep - EU
# 2. 中层蒸发计算
if WL >= params.c * params.WLM:
# 中层水量充足
EL = min(remaining_ep, WL)
new_WL = WL - EL
remaining_ep -= EL
else:
# 中层水量不足
EL = min(params.c * remaining_ep, WL)
new_WL = WL - EL
remaining_ep -= EL
# 3. 深层蒸发计算
ED = min(params.c * remaining_ep, WD)
new_WD = WD - ED
remaining_ep -= ED
# 计算总蒸发量
total_E = EU + EL + ED
# 确保含水量非负
new_WU = max(new_WU, 0)
new_WL = max(new_WL, 0)
new_WD = max(new_WD, 0)
return total_E, EU, EL, ED, new_WU, new_WL, new_WD
def calculate_runoff(PE, WU, WL, WD, params):
"""
蓄满产流计算
返回: (产流量, 上层新含水量, 中层新含水量, 深层新含水量)
"""
# 当前土壤总含水量
W = WU + WL + WD
# 无有效降雨时返回零产流
if PE <= 0:
return 0.0, WU, WL, WD
# 完全饱和情况
if W >= params.WM:
return PE, WU, WL, WD
# 计算张力水蓄水容量曲线参数
WMM = params.WM * (1 + params.b) # 最大蓄水容量
ratio = min(W / params.WM, 1.0) # 防止超过100%
A = WMM * (1 - (1 - ratio) ** (1 / (1 + params.b)))
# 产流量计算
if PE + A < WMM:
term = max(0, 1 - (PE + A) / WMM) # 防止负值
R = PE + W - params.WM + params.WM * (term ** (1 + params.b))
else:
R = PE + W - params.WM
# 确保产流在合理范围内
R = max(min(R, PE), 0)
# 计算净入渗量
net_infiltration = PE - R
# 更新土壤含水量 (分层补充)
# 1. 上层补充
add_WU = min(params.WUM - WU, net_infiltration)
new_WU = WU + add_WU
net_infiltration -= add_WU
# 2. 中层补充
add_WL = min(params.WLM - WL, net_infiltration)
new_WL = WL + add_WL
net_infiltration -= add_WL
# 3. 深层补充
add_WD = min(params.WDM - WD, net_infiltration)
new_WD = WD + add_WD
return R, new_WU, new_WL, new_WD
# ====================== 数据准备 ======================
def load_data(file_path):
"""加载输入数据"""
data = pd.read_excel(file_path)
P = data['P'].values
E0 = data['E0'].values
return P, E0
# ====================== 模拟计算主函数 ======================
def run_hydrological_model(params, P, E0):
"""运行水文模型"""
# 预处理:划分不透水面积和透水面积
P_im = P * params.Im
P_soil = P * (1 - params.Im)
# 初始化状态变量
WU, WL, WD = params.init_WU, params.init_WL, params.init_WD
# 存储结果
results = []
total_runoff = 0
total_evap = 0
for i in range(len(P)):
# 当前时段数据
current_P_total = P[i]
current_E0 = E0[i]
current_P_im = P_im[i]
current_P_soil = P_soil[i]
try:
# 三层蒸散发计算
total_E, EU, EL, ED, new_WU, new_WL, new_WD = calculate_evapotranspiration(
current_P_soil, current_E0, WU, WL, WD, params
)
# 计算净雨量
PE = current_P_soil - total_E
# 产流计算
R_soil, updated_WU, updated_WL, updated_WD = calculate_runoff(
PE, new_WU, new_WL, new_WD, params
)
# 总产流 = 透水面积产流 + 不透水面积直接径流
R_total = R_soil + current_P_im
total_runoff += R_total
total_evap += total_E
# 记录结果
result = {
'Period': i + 1,
'P_total': current_P_total,
'P_im': current_P_im,
'P_soil': current_P_soil,
'E0': current_E0,
'Ep': current_E0 * params.Kc,
'EU': EU,
'EL': EL,
'ED': ED,
'E': total_E,
'PE': PE,
'R_soil': R_soil,
'R_im': current_P_im,
'R_total': R_total,
'WU': updated_WU,
'WL': updated_WL,
'WD': updated_WD,
'Total_W': updated_WU + updated_WL + updated_WD
}
results.append(result)
# 更新状态变量
WU, WL, WD = updated_WU, updated_WL, updated_WD
except Exception as e:
print(f"错误发生在第 {i+1} 时段: {str(e)}")
print(f"当前参数: WU={WU}, WL={WL}, WD={WD}, P={current_P_total}, E0={current_E0}")
break
# 转换为DataFrame
results_df = pd.DataFrame(results)
# 添加累计量
results_df['累计产流'] = results_df['R_total'].cumsum()
results_df['累计蒸发'] = results_df['E'].cumsum()
return results_df, total_runoff, total_evap
# ====================== 结果分析与可视化 ======================
def analyze_and_visualize(results_df, params, output_dir):
"""分析结果并生成可视化图表"""
# 创建输出目录
os.makedirs(output_dir, exist_ok=True)
# 保存结果到Excel
output_path = os.path.join(output_dir, '水文模型模拟结果.xlsx')
results_df.to_excel(output_path, index=False)
print(f"模拟结果已保存到: {output_path}")
# 水量平衡验证
initial_water = params.init_WU + params.init_WL + params.init_WD
final_water = results_df.iloc[-1]['Total_W']
total_rain = results_df['P_total'].sum()
total_evap = results_df['E'].sum()
total_runoff = results_df['R_total'].sum()
balance = total_rain - (final_water - initial_water) - total_evap - total_runoff
print("\n水量平衡验证:")
print(f"总降雨量: {total_rain:.2f} mm")
print(f"总蒸发量: {total_evap:.2f} mm")
print(f"总径流量: {total_runoff:.2f} mm")
print(f"土壤蓄水变化: {final_water - initial_water:.2f} mm")
print(f"水量平衡误差: {balance:.6f} mm (应接近0)")
# 可视化
plt.figure(figsize=(15, 10))
# 1. 降雨-径流过程
plt.subplot(2, 2, 1)
plt.plot(results_df['Period'], results_df['P_total'], 'b-', label='降雨量')
plt.plot(results_df['Period'], results_df['R_total'], 'r-', label='径流量')
plt.xlabel('时段')
plt.ylabel('水量 (mm)')
plt.title('降雨-径流过程线')
plt.legend()
plt.grid(True)
# 2. 土壤含水量变化
plt.subplot(2, 2, 2)
plt.plot(results_df['Period'], results_df['WU'], 'g-', label='上层土壤水')
plt.plot(results_df['Period'], results_df['WL'], 'b-', label='中层土壤水')
plt.plot(results_df['Period'], results_df['WD'], 'r-', label='深层土壤水')
plt.xlabel('时段')
plt.ylabel('土壤含水量 (mm)')
plt.title('土壤含水量变化过程')
plt.legend()
plt.grid(True)
# 3. 累计水量变化
plt.subplot(2, 2, 3)
plt.plot(results_df['Period'], results_df['累计产流'], 'r-', label='累计径流')
plt.plot(results_df['Period'], results_df['累计蒸发'], 'b-', label='累计蒸发')
plt.xlabel('时段')
plt.ylabel('累计水量 (mm)')
plt.title('累计水量变化过程')
plt.legend()
plt.grid(True)
# 4. 蒸发分量构成
plt.subplot(2, 2, 4)
plt.stackplot(
results_df['Period'],
results_df['EU'],
results_df['EL'],
results_df['ED'],
labels=['上层蒸发', '中层蒸发', '深层蒸发']
)
plt.xlabel('时段')
plt.ylabel('蒸发量 (mm)')
plt.title('蒸发分量构成')
plt.legend(loc='upper left')
plt.grid(True)
plt.tight_layout()
# 保存图表
plot_path = os.path.join(output_dir, '水文模型结果可视化.png')
plt.savefig(plot_path)
print(f"可视化结果已保存到: {plot_path}")
plt.close()
# ====================== 主程序 ======================
if __name__ == "__main__":
# 记录开始时间
start_time = time.time()
# 初始化模型参数
params = ModelParams()
# 加载数据
data_file = r"C:\Users\80401\Desktop\SHUJU2.xlsx"
P, E0 = load_data(data_file)
# 运行模型
results_df, total_runoff, total_evap = run_hydrological_model(params, P, E0)
# 分析结果并可视化
output_dir = os.path.join(os.path.expanduser('~'), 'Desktop', '水文模型结果')
analyze_and_visualize(results_df, params, output_dir)
# 计算执行时间
end_time = time.time()
print(f"\n模型执行完成,耗时: {end_time - start_time:.2f}秒")
print(f"总径流量: {total_runoff:.2f} mm")
print(f"总蒸发量: {total_evap:.2f} mm")
修正下列错误:D:\Python\python.exe D:\Python\PythonProject\EPW\a.py
Traceback (most recent call last):
File "D:\Python\PythonProject\EPW\a.py", line 122, in <module>
R_soil, updated_WU, updated_WL, updated_WD = calculate_runoff(
~~~~~~~~~~~~~~~~^
PE, new_WU, new_WL, new_WD
^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "D:\Python\PythonProject\EPW\a.py", line 55, in calculate_runoff
if PE <= 0:
^^^^^^^
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
进程已结束,退出代码为 1
最新发布