1. 项目概述用Python解放水文地质计算干了这么多年水文地质最头疼的就是每次做完抽水试验抱着一堆现场记录数据回来在Excel里吭哧吭哧套公式、查表、画图一个参数算错后面全得重来。特别是计算渗透系数K和影响半径R这种核心参数公式复杂迭代计算多手动处理不仅效率低还容易出错。后来我开始用Python把这些计算流程自动化才发现这才是“解放生产力”的正确姿势。这个项目就是把我这些年用Python处理承压水完整井稳定流抽水试验数据、计算K和R的经验整理成一套清晰、可复现的代码流程。它不是什么高深的AI模型而是一个解决具体工程问题的实用工具箱。核心目标是输入现场观测的降深、流量、时间等原始数据自动完成公式套用、迭代求解、结果可视化并输出规范的计算报告。无论是刚入行的技术员还是想提升效率的工程师都能通过这套代码把繁琐的计算工作交给计算机自己则专注于更重要的数据分析与地质判断。2. 核心原理与公式拆解理解计算背后的地质逻辑在动手写代码之前我们必须吃透背后的水文地质原理。这不是简单的数学计算每一个参数都对应着含水层的物理特性。2.1 承压水完整井稳定流抽水的模型假设我们讨论的“承压水完整井稳定流抽水”是基于泰斯Theis公式的简化版本——裘布依Dupuit公式。它有几个关键假设理解这些是正确应用的前提含水层均质、等厚、水平无限延伸这是理想情况实际含水层总有差异但计算时我们以此为基础。抽水前地下水处于稳定状态初始水头面是水平的。水流服从达西定律即渗流速度与水力坡度成正比。井为完整井井的滤管贯穿整个含水层厚度。抽水流量恒定从开始到结束抽水机的出水量Q保持不变。水流为稳定流抽水一段时间后观测井中的水位降深s不再随时间变化达到稳定状态。注意现场中“稳定”是相对的通常指单位时间内的水位变幅小于某个阈值如5cm/h。我们的计算依赖于这个“稳定”时刻的数据。2.2 渗透系数K与影响半径R的计算公式基于上述模型有两个核心公式1. 渗透系数K的计算公式裘布依公式这是最核心的公式用于计算含水层导水能力。K (Q / (2πM s)) * ln(R/r)K: 渗透系数m/d这是我们要求解的核心参数之一。Q: 抽水井的稳定流量m³/d。M: 承压含水层的厚度m。s: 在距离抽水井r处的观测井或抽水井自身的稳定水位降深m。R: 影响半径m即抽水影响范围的半径。r: 观测井到抽水井的距离m。如果是用抽水井自身的降深则r为抽水井的半径rw。ln: 自然对数。2. 影响半径R的经验公式库萨金公式影响半径R在公式中与K耦合无法直接分离求解。实践中常用经验公式进行估算其中库萨金公式较为常用R 10 * s * sqrt(K)s: 抽水井的水位降深m。K: 渗透系数m/d。看这里又需要K所以K和R是互相依赖的。3. 迭代求解的逻辑闭环从上面两个公式可以看出我们陷入了一个“先有鸡还是先有蛋”的循环算K需要R算R需要K。因此必须采用迭代法求解。首先给影响半径R一个初始估计值R0例如根据经验取100m或200m。将R0代入裘布依公式计算出第一个渗透系数K1。将计算出的K1代入库萨金公式得到一个新的影响半径R1。比较R1与R0的差值。如果差值大于我们设定的容差如0.01m则将R1作为新的R0重复步骤2-3。如此反复直到两次计算出的R值之差小于容差此时对应的K和R即为最终解。这个过程手动计算极其繁琐但正是Python等编程语言的用武之地。3. 开发环境搭建与核心工具选型工欲善其事必先利其器。一个高效、清晰的开发环境能让后续编码事半功倍。3.1 Python环境与IDE选择我强烈推荐使用Anaconda来管理Python环境。它集成了科学计算所需的绝大部分库如NumPy, Pandas, Matplotlib并且环境隔离做得很好避免项目间的包版本冲突。安装Anaconda从官网下载安装包一路默认安装即可。安装后你会拥有一个基础的base环境。为项目创建独立环境打开Anaconda PromptWindows或终端Mac/Linux执行以下命令创建一个名为hydro_calc的新环境并指定Python版本如3.9conda create -n hydro_calc python3.9 conda activate hydro_calcIDE选择VS Code是目前最灵活轻量的选择。安装Python扩展后智能提示、调试、代码管理功能都非常强大。当然如果你习惯了PyCharm的专业版那也是极好的选择。关键在于顺手。3.2 必需的三方库及其作用在我们的环境中需要安装以下几个核心库# 在激活的 hydro_calc 环境中执行 pip install numpy pandas matplotlib scipyNumPy进行高效的数组运算和数学计算比如对数运算、迭代循环。Pandas数据处理的核心。用于读取、清洗、组织我们的抽水试验数据如多个时间点、多个观测井的降深数据功能比Excel强大得多。Matplotlib绘制专业图表如降深-时间曲线s-t曲线、降深-距离曲线s-r曲线可视化是验证数据稳定性和展示成果的关键。SciPy高级科学计算库。虽然本次基础计算用不到其复杂功能但其优化模块在未来处理非稳定流或复杂模型时会很有用。实操心得务必在项目开始时就用conda或pip生成一个requirements.txt文件pip freeze requirements.txt记录所有库的精确版本。这样在任何其他电脑上重建环境时只需pip install -r requirements.txt即可完美复现避免“在我电脑上好好的”这种问题。4. 数据准备与预处理从野外记录到规整数据框原始数据通常来自野外记录本或Excel格式可能杂乱。Python计算的第一步就是将其转化为程序可读的、规整的结构。4.1 设计合理的数据结构我通常建议用一个字典或Pandas的DataFrame来组织一次抽水试验的所有参数和数据。以下是一个示例数据结构# 示例定义抽水试验的基本参数和观测数据 test_data { # 1. 基本水文地质参数 aquifer_thickness: 25.0, # M含水层厚度 (m) pumping_rate: 1200.0, # Q抽水流量 (m³/d) well_radius: 0.2, # rw抽水井半径 (m) # 2. 观测井信息列表支持多个观测井 observation_wells: [ {name: OB-1, distance: 10.0, stable_drawdown: 3.2}, # r (m), s (m) {name: OB-2, distance: 25.0, stable_drawdown: 1.8}, {name: OB-3, distance: 50.0, stable_drawdown: 0.9}, ], # 3. 抽水井自身降深可选用于计算R的初始估算 main_well_drawdown: 5.5, # s_w (m) # 4. 迭代计算参数 initial_R_guess: 200.0, # R的初始猜测值 (m) tolerance: 0.01, # 迭代收敛容差 (m) max_iterations: 100 # 最大迭代次数防止无限循环 }4.2 数据读取与清洗数据可能来自CSV文件。我们需要用Pandas读取并检查。import pandas as pd # 假设有一个‘drawdown_data.csv’文件列包括Time(min), OB1_s(m), OB2_s(m), OB3_s(m) df pd.read_csv(drawdown_data.csv) # 查看数据前几行和基本信息 print(df.head()) print(df.info()) # 数据清洗检查缺失值 if df.isnull().sum().any(): print(发现缺失值需要处理。) # 可以根据前后数据插值或删除该时间点谨慎 # df df.interpolate() # 线性插值 # df df.dropna() # 删除含缺失值的行 # 关键步骤确定稳定降深s # 通常我们取时间序列后期波动较小的平均值作为稳定降深 # 例如取最后30分钟的数据 stable_period df[df[Time(min)] df[Time(min)].max() - 30] stable_drawdowns stable_period[[OB1_s(m), OB2_s(m), OB3_s(m)]].mean() print(各观测井的稳定降深m:) print(stable_drawdowns)注意事项确定“稳定降深”是计算的关键也是容易出错的地方。不能简单地取最后一个值。一定要绘制s-t曲线肉眼判断水位是否进入平稳阶段并取该阶段的平均值。自动化脚本可以设定一个阈值如连续10分钟降深变化率小于0.5%但人工复核必不可少。5. 核心算法实现迭代求解K与R有了干净的数据我们就可以实现第2章中提到的迭代算法了。我们将这个过程封装成一个函数提高代码的复用性和可读性。5.1 单观测井迭代计算函数首先我们实现针对一个观测井数据计算K和R的函数。import numpy as np def calculate_k_r_for_one_well(Q, M, s, r, s_mainNone, R_initial200.0, tol0.01, max_iter100): 根据裘布依公式和库萨金公式迭代计算渗透系数K和影响半径R。 参数: Q: 抽水流量 (m³/d) M: 含水层厚度 (m) s: 观测井稳定降深 (m) r: 观测井到抽水井距离 (m) s_main: 抽水井自身降深 (m)用于库萨金公式。如果为None则使用观测井降深s估算。 R_initial: 影响半径初始猜测值 (m) tol: 迭代收敛容差 (m) max_iter: 最大迭代次数 返回: K, R, iterations, converged K: 渗透系数 (m/d) R: 影响半径 (m) iterations: 实际迭代次数 converged: 是否收敛 (布尔值) R_old R_initial if s_main is None: s_main s # 若无主井降深则用观测井降深近似 for i in range(max_iter): # 1. 使用当前的R_old通过裘布依公式计算K # 公式: K (Q / (2 * π * M * s)) * ln(R/r) # 注意对数内R/r必须大于1否则无物理意义 if R_old r: raise ValueError(f迭代过程中R({R_old}) r({r})无物理意义请检查初始值或数据。) K_new (Q / (2 * np.pi * M * s)) * np.log(R_old / r) # 2. 使用计算出的K_new通过库萨金公式更新R # 公式: R 10 * s * sqrt(K) R_new 10 * s_main * np.sqrt(K_new) # 3. 检查是否收敛 if abs(R_new - R_old) tol: return K_new, R_new, i1, True # 4. 未收敛更新R_old继续迭代 R_old R_new # 如果达到最大迭代次数仍未收敛 print(f警告未在{max_iter}次迭代内收敛。最后计算的K{K_new:.4f}, R{R_new:.4f}) return K_new, R_new, max_iter, False5.2 多观测井数据处理与结果整合一次抽水试验通常有多个观测井每个井都可以算出一组K、R。理论上在理想条件下它们应该相近。我们可以计算多组结果并求平均值或分析其离散程度这本身就是对试验质量的一种检验。def calculate_from_multiple_wells(test_data): 处理包含多个观测井的试验数据并汇总结果。 Q test_data[pumping_rate] M test_data[aquifer_thickness] s_main test_data.get(main_well_drawdown) # 获取主井降深可能为None results [] for well in test_data[observation_wells]: name well[name] r well[distance] s well[stable_drawdown] try: K, R, iter_cnt, converged calculate_k_r_for_one_well( Q, M, s, r, s_main, R_initialtest_data[initial_R_guess], toltest_data[tolerance], max_itertest_data[max_iterations] ) results.append({ 观测井: name, 距离r(m): r, 降深s(m): s, 渗透系数K(m/d): K, 影响半径R(m): R, 迭代次数: iter_cnt, 是否收敛: 是 if converged else 否 }) print(f{name}: K {K:.4f} m/d, R {R:.2f} m (迭代{iter_cnt}次)) except ValueError as e: print(f计算{name}时出错{e}) results.append({ 观测井: name, 距离r(m): r, 降深s(m): s, 渗透系数K(m/d): None, 影响半径R(m): None, 迭代次数: None, 是否收敛: 错误 }) # 将结果转换为DataFrame便于分析 results_df pd.DataFrame(results) print(\n 计算结果汇总 ) print(results_df.to_string(indexFalse)) # 计算平均K和R排除计算错误的井 valid_results results_df[results_df[是否收敛].isin([是])] if not valid_results.empty: avg_K valid_results[渗透系数K(m/d)].mean() avg_R valid_results[影响半径R(m)].mean() std_K valid_results[渗透系数K(m/d)].std() print(f\n基于{len(valid_results)}个有效观测井) print(f平均渗透系数 K_avg {avg_K:.4f} ± {std_K:.4f} m/d) print(f平均影响半径 R_avg {avg_R:.2f} m) # 可以进一步分析例如离散系数标准差/平均值是否过大判断含水层均质性 if std_K / avg_K 0.2: # 假设离散系数大于20%提示不均质 print( 注意各观测井计算的K值离散较大含水层可能非均质或试验数据存在问题。) return results_df运行这个函数输入之前准备好的test_data字典就能得到一份清晰的成果汇总表。6. 结果可视化与报告生成数字结果虽然精确但图表更能直观揭示规律和问题。同时生成一份格式规范的计算报告是交付成果的必要步骤。6.1 绘制专业图表至少需要绘制两种核心图件1. 降深-距离关系图s-r图在双对数坐标纸上稳定流理论的s与ln(r)应是线性关系。绘制此图可以直观验证数据是否符合理论模型。import matplotlib.pyplot as plt def plot_drawdown_distance(results_df, test_data): 绘制降深s与距离r的关系图半对数坐标s-ln(r)。 fig, ax plt.subplots(figsize(8, 6)) r_vals results_df[距离r(m)] s_vals results_df[降深s(m)] # 绘制散点 ax.semilogx(r_vals, s_vals, bo-, linewidth1.5, markersize8, label观测数据) # 根据裘布依公式s与ln(r)呈线性关系s (Q/(2πKM)) * ln(R/r) # 我们可以用计算出的平均K和R来绘制理论曲线 Q test_data[pumping_rate] M test_data[aquifer_thickness] # 假设我们取平均K和R avg_K results_df[渗透系数K(m/d)].mean() avg_R results_df[影响半径R(m)].mean() # 生成一系列r值用于绘制平滑曲线 r_smooth np.logspace(np.log10(min(r_vals)*0.5), np.log10(avg_R*1.2), 100) s_theory (Q / (2 * np.pi * avg_K * M)) * np.log(avg_R / r_smooth) ax.semilogx(r_smooth, s_theory, r--, linewidth2, labelf理论曲线 (K{avg_K:.2f}, R{avg_R:.0f})) ax.set_xlabel(距离抽水井距离 r (m) - 对数坐标, fontsize12) ax.set_ylabel(稳定降深 s (m), fontsize12) ax.set_title(承压完整井稳定流抽水试验 s-ln(r) 关系图, fontsize14) ax.grid(True, whichboth, linestyle--, alpha0.6) ax.legend() plt.tight_layout() plt.savefig(s_r_relationship.png, dpi300) # 保存高清图 plt.show()2. 各观测井计算参数对比图用柱状图或散点图展示各井计算出的K、R值一目了然地看出其一致性和离散度。def plot_parameter_comparison(results_df): 绘制各观测井计算出的K和R值对比图。 fig, (ax1, ax2) plt.subplots(1, 2, figsize(12, 5)) wells results_df[观测井] K_vals results_df[渗透系数K(m/d)] R_vals results_df[影响半径R(m)] # 渗透系数K对比 bars1 ax1.bar(wells, K_vals, colorskyblue, edgecolorblack) ax1.axhline(yK_vals.mean(), colorred, linestyle--, labelf平均值: {K_vals.mean():.3f}) ax1.set_ylabel(渗透系数 K (m/d), fontsize12) ax1.set_title(各观测井计算渗透系数对比, fontsize14) ax1.legend() # 在柱子上方标注数值 for bar in bars1: height bar.get_height() ax1.text(bar.get_x() bar.get_width()/2., height 0.001, f{height:.3f}, hacenter, vabottom, fontsize9) # 影响半径R对比 bars2 ax2.bar(wells, R_vals, colorlightgreen, edgecolorblack) ax2.axhline(yR_vals.mean(), colorred, linestyle--, labelf平均值: {R_vals.mean():.1f}) ax2.set_ylabel(影响半径 R (m), fontsize12) ax2.set_title(各观测井计算影响半径对比, fontsize14) ax2.legend() for bar in bars2: height bar.get_height() ax2.text(bar.get_x() bar.get_width()/2., height 1, f{height:.0f}, hacenter, vabottom, fontsize9) plt.tight_layout() plt.savefig(parameter_comparison.png, dpi300) plt.show()6.2 生成计算报告最后我们可以将输入参数、计算过程、结果和图表整合成一份文本报告。def generate_report(test_data, results_df, avg_K, avg_R): 生成简单的文本计算报告。 report_lines [] report_lines.append(*60) report_lines.append( 承压水完整井稳定流抽水试验计算报告) report_lines.append(*60) report_lines.append(f\n一、试验基本参数) report_lines.append(f 抽水流量 Q {test_data[pumping_rate]} m³/d) report_lines.append(f 含水层厚度 M {test_data[aquifer_thickness]} m) report_lines.append(f 抽水井半径 rw {test_data[well_radius]} m) report_lines.append(f 抽水井降深 sw {test_data.get(main_well_drawdown, 未提供)} m) report_lines.append(f\n二、观测井数据) for well in test_data[observation_wells]: report_lines.append(f {well[name]}: 距离 r {well[distance]} m, 稳定降深 s {well[stable_drawdown]} m) report_lines.append(f\n三、迭代计算参数) report_lines.append(f 初始影响半径猜测 R0 {test_data[initial_R_guess]} m) report_lines.append(f 收敛容差 tolerance {test_data[tolerance]} m) report_lines.append(f 最大迭代次数 max_iterations {test_data[max_iterations]}) report_lines.append(f\n四、各观测井计算结果) report_lines.append(results_df.to_string(indexFalse)) report_lines.append(f\n五、推荐参数值基于有效观测井平均值) report_lines.append(f 建议采用的渗透系数 K {avg_K:.4f} m/d) report_lines.append(f 建议采用的影响半径 R {avg_R:.2f} m) report_lines.append(f\n六、备注) report_lines.append( 1. 计算基于裘布依稳定流公式及库萨金经验公式。) report_lines.append( 2. 结果适用于均质、等厚、无限延伸的理想承压含水层。) report_lines.append( 3. 实际工程应用需结合地质条件进行综合判断。) report_lines.append(\n *60) report_lines.append(报告生成完成。) report_text \n.join(report_lines) # 保存报告到文件 with open(pumping_test_calculation_report.txt, w, encodingutf-8) as f: f.write(report_text) print(report_text) # 同时在控制台输出 return report_text7. 常见问题、误差分析与实战心得在实际应用中你一定会遇到各种问题。下面是我踩过的一些坑和对应的解决思路。7.1 迭代计算不收敛或结果异常问题现象程序报错“R r”或迭代几百次也不收敛或计算出的K值极大/极小。原因排查数据输入错误检查Q、M、s、r的单位是否统一强烈建议全部转换为“米-天”制。流量Q是m³/d不是m³/h降深s是米不是厘米。初始R猜测值不合理R_initial不能小于观测井距离r。如果观测井很远如200m初始猜测值至少要比它大。可以尝试根据经验公式R ≈ 3000 * s * sqrt(K)更粗略先估算一个数量级或者直接设一个较大的值如500m、1000m。降深s为0或极小如果降深测量误差导致s接近0公式中分母趋近于0K会趋于无穷大。需要检查观测数据是否可靠。含水层非均质性强烈实际条件严重偏离“均质”假设导致公式本身不适用。此时不同观测井计算的K值会差异巨大。解决方案在calculate_k_r_for_one_well函数中增加更严格的输入校验。添加一个“安全模式”当迭代超过一定次数或R值异常波动时自动调整初始猜测值或终止计算并给出明确警告。绘制s-ln(r)图。如果数据点明显偏离直线则提示用户理论模型可能不适用。7.2 多个观测井计算结果离散度大问题现象三个观测井算出的K分别是1.2, 5.6, 0.8 m/d相差数倍。地质含义这很可能揭示了含水层的非均质性。距离近的井可能受到局部裂隙或夹层的影响。处理建议不要简单取平均应分析离散原因。是某个井的数据异常如s未真正稳定还是地质条件确实如此分区给出参数如果含水层有明显分区如上下游可以分区统计K值。在报告中明确说明给出平均值的同时必须注明标准差和离散系数并附上“计算结果离散度较大建议结合地质勘察资料综合分析”的说明。这是专业性的体现。7.3 如何选择最终推荐值这是工程判断而不仅仅是数学计算。优先考虑距离抽水井适中的观测井太近的井可能受井损影响太远的井降深小、测量相对误差大。通常认为1.5倍含水层厚度以外的观测井数据更可靠。参考抽水井自身降深计算的结果如果抽水井的降深数据质量高用rrw和ssw计算出的K值具有重要参考意义。与地区经验值对比将计算结果与同一地区、同类地层的经验渗透系数范围进行对比如果偏离太远需要回头检查数据。保守原则对于涉及安全的设计如基坑降水在参数离散时有时会倾向于选取偏不利如较大的K值进行设计。7.4 代码优化与扩展方向当这个基础工具用顺手后你可以考虑以下扩展图形用户界面GUI使用PyQt或Tkinter打包成一个桌面小软件方便野外技术人员直接输入Excel数据点按钮出结果。非稳定流计算集成泰斯Theis公式或雅各布Jacob近似公式处理更普通的非稳定流抽水试验数据这需要用到scipy.special中的指数积分函数。自动识别稳定段编写算法自动从s-t时间序列数据中识别出水位稳定阶段并提取s值实现全流程自动化。生成Word/PDF报告使用python-docx或ReportLab库将文字、表格、图片自动排版生成可直接交付的正式报告文档。这套代码的终极价值在于它将你从重复、易错的手工计算中解放出来让你有更多时间去思考数据背后的地质故事。一开始搭建框架会花点时间但一旦建成它就是你的专属“数字助手”所有同类项目的计算效率都能提升十倍以上。