import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
from scipy.stats import norm
import tkinter as tk
from tkinter import filedialog
import os
def read_section_file(file_path):
"""读取断面数据文件"""
data = np.loadtxt(file_path)
s = data[:, 0] # 起点距
z = data[:, 1] # 高程
n = 0.03 # 默认糙率
return s, z, n
def project_section(s, z, angle):
"""断面投影变换（考虑桥梁斜交角度）"""
rad = np.deg2rad(angle)
s_proj = s * np.cos(rad)
return s_proj, z
def calc_area_perimeter(s, z, water_level):
"""计算过水断面面积和湿周"""
# 找到水位与断面的交点
s_under = []
z_under = []
# 找出所有低于水位的点
for i in range(len(s) - 1):
if z[i] <= water_level or z[i + 1] <= water_level:
# 左交点
if z[i] > water_level and z[i + 1] <= water_level:
ratio = (water_level - z[i]) / (z[i + 1] - z[i])
s_left = s[i] + ratio * (s[i + 1] - s[i])
s_under.append(s_left)
z_under.append(water_level)
# 右交点
elif z[i] <= water_level and z[i + 1] > water_level:
ratio = (water_level - z[i]) / (z[i + 1] - z[i])
s_right = s[i] + ratio * (s[i + 1] - s[i])
s_under.append(s_right)
z_under.append(water_level)
# 完全在水下的点
if z[i] <= water_level:
s_under.append(s[i])
z_under.append(z[i])
# 如果没有交点，返回0
if len(s_under) < 2:
return 0, 0
# 按起点距排序
idx = np.argsort(s_under)
s_under = np.array(s_under)[idx]
z_under = np.array(z_under)[idx]
# 计算过水面积（梯形法）
area = 0
for i in range(len(s_under) - 1):
width = s_under[i + 1] - s_under[i]
depth_left = water_level - z_under[i]
depth_right = water_level - z_under[i + 1]
area += width * (depth_left + depth_right) / 2
# 计算湿周（河床长度）
perimeter = 0
for i in range(len(s_under) - 1):
dx = s_under[i + 1] - s_under[i]
dy = z_under[i + 1] - z_under[i]
perimeter += np.sqrt(dx ** 2 + dy ** 2)
return area, perimeter
def get_G(s, z, water_level, n=0.03):
"""计算流量模数G"""
area, perimeter = calc_area_perimeter(s, z, water_level)
if area <= 0 or perimeter <= 0:
return 0
R = area / perimeter # 水力半径
G = area * R ** (2 / 3) / n # 曼宁公式
return G
def get_Kp(Cv, Cs2Cv, freq):
"""计算P-III型分布的模比系数Kp"""
# 使用正态分布近似
p = 1 - freq
Kp = 1 + norm.ppf(p) * Cv
return Kp
def save_result(file_path, design_water, safe_water, warn_water,
K0, K1, K2, G_design, G_safe, G_warn, z_curve, G_curve):
"""保存计算结果到文件"""
with open(file_path, 'w') as f:
f.write("桥梁洪水位计算结果\n")
f.write("=" * 50 + "\n")
f.write(f"设计水位: {design_water:.2f} m\n")
f.write(f"安全运行水位: {safe_water:.2f} m\n")
f.write(f"警戒水位: {warn_water:.2f} m\n")
f.write(f"设计流量模数: {G_design:.2f} m³/s\n")
f.write(f"安全流量模数: {G_safe:.2f} m³/s\n")
f.write(f"警戒流量模数: {G_warn:.2f} m³/s\n")
f.write(f"模比系数K0: {K0:.4f}\n")
f.write(f"模比系数K1: {K1:.4f}\n")
f.write(f"模比系数K2: {K2:.4f}\n")
f.write("\n水位-流量模数曲线:\n")
f.write("水位(m)\t流量模数(m³/s)\n")
for z, g in zip(z_curve, G_curve):
f.write(f"{z:.2f}\t{g:.2f}\n")
def main():
"""主函数"""
# 创建Tkinter根窗口（不显示）
root = tk.Tk()
root.withdraw()
# 文件选择对话框
design_file = filedialog.askopenfilename(
title="选择设计断面数据文件",
filetypes=[("文本文件", "*.txt"), ("所有文件", "*.*")]
)
current_file = filedialog.askopenfilename(
title="选择实测断面数据文件",
filetypes=[("文本文件", "*.txt"), ("所有文件", "*.*")]
)
save_file = filedialog.asksaveasfilename(
title="保存计算结果",
filetypes=[("文本文件", "*.txt")],
defaultextension=".txt"
)
# 用户输入参数
dAngle = float(input("请输入大桥斜交角度(度): ") or "90")
Cv = float(input("请输入洪峰流量变差系数(0.2-0.6): ") or "0.35")
Cs2Cv = float(input("请输入偏态系数与变差系数之比(2-4): ") or "3")
dKnowWaterLevel = float(input("请输入设计水位(m): ") or "64")
dKnowFreq = 1 / float(input("请输入设计水位对应的频率(如100表示1/100): ") or "100")
dSafeFreq = 1 / float(input("请输入安全运行水位对应的频率: ") or "100")
dWarnFreq = 1 / float(input("请输入警戒水位对应的频率: ") or "25")
# 读取断面数据
dDs, dDz, dDn = read_section_file(design_file)
dMs, dMz, dMn = read_section_file(current_file)
# 断面投影
dDs, dDz = project_section(dDs, dDz, dAngle)
dMs, dMz = project_section(dMs, dMz, dAngle)
# 绘制断面图
plt.figure(figsize=(10, 6))
plt.plot(dDs, dDz, 'b-', linewidth=2, label='设计断面')
plt.plot(dMs, dMz, 'r-', linewidth=2, label='实测断面')
plt.xlabel('起点距(m)')
plt.ylabel('高程(m)')
plt.title('桥梁断面图')
plt.legend()
plt.grid(True)
plt.show()
# 计算已知水位对应的流量模数
G_Know = get_G(dDs, dDz, dKnowWaterLevel, dDn)
# 计算水位-流量模数曲线
minMz, maxMz = min(dMz), max(dMz)
minMs, maxMs = min(dMs), max(dMs)
U = 100 # 曲线点数
dCurveZ = np.linspace(minMz, maxMz, U)
dCurveG = np.zeros(U)
for i, z in enumerate(dCurveZ):
dCurveG[i] = get_G(dMs, dMz, z, dMn)
# 计算模比系数
dK0 = get_Kp(Cv, Cs2Cv, dKnowFreq)
dK1 = get_Kp(Cv, Cs2Cv, dSafeFreq)
dK2 = get_Kp(Cv, Cs2Cv, dWarnFreq)
# 计算安全水位和警戒水位
G_Safe = G_Know * dK1 / dK0
G_Warn = G_Know * dK2 / dK0
# 插值计算水位
f = interpolate.interp1d(dCurveG, dCurveZ, bounds_error=False, fill_value="extrapolate")
dSafeWaterLevel = f(G_Safe)
dWarnWaterLevel = f(G_Warn)
# 输出结果
print(f"设计水位 = {dKnowWaterLevel:.2f}m")
print(f"安全运行水位 = {dSafeWaterLevel:.2f}m")
print(f"警戒水位 = {dWarnWaterLevel:.2f}m")
print(f"安全水位流量模数 = {G_Safe:.2f}m³/s")
print(f"警戒水位流量模数 = {G_Warn:.2f}m³/s")
# 绘制结果图
plt.figure(figsize=(12, 8))
plt.plot(dCurveG, dCurveZ, 'k-', linewidth=2, label='水位-流量模数曲线')
plt.axvline(x=G_Safe, color='r', linestyle='-.', linewidth=2)
plt.axhline(y=dSafeWaterLevel, color='r', linestyle='-.', linewidth=2)
plt.axvline(x=G_Warn, color='b', linestyle=':', linewidth=2)
plt.axhline(y=dWarnWaterLevel, color='b', linestyle=':', linewidth=2)
plt.annotate(f'安全运行基准洪水流量模数\n$G_p$={G_Safe:.2f} m³/s',
xy=(G_Safe, dSafeWaterLevel),
xytext=(G_Safe * 1.05, dSafeWaterLevel * 0.9),
arrowprops=dict(arrowstyle='->'), fontsize=12)
plt.annotate(f'安全运行基准水位\n$h_p$={dSafeWaterLevel:.2f} m',
xy=(G_Safe, dSafeWaterLevel),
xytext=(G_Safe * 0.5, dSafeWaterLevel * 1.05),
fontsize=12)
plt.annotate(f'警戒洪水流量模数\n$G_y$={G_Warn:.2f} m³/s',
xy=(G_Warn, dWarnWaterLevel),
xytext=(G_Warn * 0.95, dWarnWaterLevel * 0.85),
arrowprops=dict(arrowstyle='->'), fontsize=12,
horizontalalignment='right')
plt.annotate(f'警戒水位\n$h_y$={dWarnWaterLevel:.2f} m',
xy=(G_Warn, dWarnWaterLevel),
xytext=(G_Warn * 0.5, dWarnWaterLevel * 0.95),
fontsize=12)
plt.xlabel('流量模数 $G$ (m³/s)', fontsize=14, fontweight='bold')
plt.ylabel('水位 $h$ (m)', fontsize=14, fontweight='bold')
plt.title('桥梁洪水位计算结果', fontsize=16)
plt.grid(True)
plt.xlim(min(dCurveG) * 0.9, max(dCurveG) * 1.1)
plt.ylim(min(dCurveZ) * 0.95, max(dCurveZ) * 1.05)
plt.tight_layout()
plt.show()
# 保存结果
save_result(save_file, dKnowWaterLevel, dSafeWaterLevel, dWarnWaterLevel,
dK0, dK1, dK2, G_Know, G_Safe, G_Warn, dCurveZ, dCurveG)
print(f"计算结果已保存到: {save_file}")
if __name__ == "__main__":
main()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 39640 (\N{CJK UNIFIED IDEOGRAPH-9AD8}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 31243 (\N{CJK UNIFIED IDEOGRAPH-7A0B}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 26725 (\N{CJK UNIFIED IDEOGRAPH-6865}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 26753 (\N{CJK UNIFIED IDEOGRAPH-6881}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 26029 (\N{CJK UNIFIED IDEOGRAPH-65AD}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 38754 (\N{CJK UNIFIED IDEOGRAPH-9762}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 22270 (\N{CJK UNIFIED IDEOGRAPH-56FE}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 36215 (\N{CJK UNIFIED IDEOGRAPH-8D77}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 28857 (\N{CJK UNIFIED IDEOGRAPH-70B9}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 36317 (\N{CJK UNIFIED IDEOGRAPH-8DDD}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 35774 (\N{CJK UNIFIED IDEOGRAPH-8BBE}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 35745 (\N{CJK UNIFIED IDEOGRAPH-8BA1}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 23454 (\N{CJK UNIFIED IDEOGRAPH-5B9E}) missing from font(s) DejaVu Sans.
plt.show()
D:\PythonProject1\水位计算软件.py:168: UserWarning: Glyph 27979 (\N{CJK UNIFIED IDEOGRAPH-6D4B}) missing from font(s) DejaVu Sans.
plt.show()
最新发布
博客给出了一个用于计算两个给定时间差的函数。该函数接收两个时间参数，通过一系列计算得出时间差，并以“X天X小时X分X秒”的格式输出结果，方便在程序中调用。
扫一扫
1431
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
