基于MATLAB的有杆抽油系统动力学建模与智能故障诊断实践

1. 项目缘起：从“黑箱”到“白箱”的抽油系统认知跃迁

在石油开采的现场，有杆抽油系统（俗称“磕头机”）是陆地油田最常见的一道风景。这套机械系统看似结构简单，但其内部动力学行为却异常复杂。在我早期参与油田数字化改造项目时，面对一个反复出现“光杆断脱”故障的井，传统的经验诊断方法——听声音、看电流图、凭老师傅的手感——显得力不从心。我们耗费了大量时间排查，更换了悬绳器、调整了冲程冲次，问题却依旧间歇性出现，直接导致了近一个月的非计划性停产和可观的产量损失。

这次经历让我深刻意识到，仅仅依靠外部现象和经验去理解这套由地面驱动设备、抽油杆柱、井下泵及油管组成的“长链条”系统，如同隔靴搔痒。我们必须建立一套能够精确描述其内在物理规律的数学模型，将“黑箱”变为“白箱”。而MATLAB，以其强大的数值计算、矩阵处理、控制系统仿真及丰富的可视化工具箱，成为了实现这一目标的首选利器。本次分享的“有杆抽油系统的数学建模及诊断”，正是基于这样的工程背景，旨在通过严谨的数学推导和仿真实践，构建一套可用于系统性能分析、故障预测与智能诊断的数字化工具链。无论你是从事油气田开发的研究人员、现场工程师，还是对机电系统建模与控制感兴趣的学生，这套方法都能为你提供一个从理论到实践的完整视角。

2. 有杆抽油系统动力学模型的核心：一维波动方程及其离散化

有杆抽油系统的核心物理过程是能量通过细长的抽油杆柱（长度可达数千米）从地面传递到井下泵。抽油杆柱在交变载荷作用下产生的纵向振动是分析一切问题（如应力、位移、载荷）的基础。描述这一振动的最经典模型是一维波动方程。

2.1 一维波动方程的推导与物理意义

我们首先将抽油杆柱视为一个均质、连续的弹性杆。根据牛顿第二定律和胡克定律，可以推导出描述杆柱纵向振动的一维波动方程：

∂²u(x,t)/∂t² = a² * ∂²u(x,t)/∂x² - c * ∂u(x,t)/∂t

其中：

• u(x, t)是距离井口x处、时间t时杆柱的位移（米）。
• a = sqrt(E/ρ)是应力波在杆柱中的传播速度（米/秒），E是杆材的弹性模量（帕斯卡），ρ是杆材密度（千克/立方米）。对于钢杆，a通常在5000 m/s左右。
• c是粘滞阻尼系数（1/秒），用于表征杆柱在油液中的运动阻尼。

注意：这个方程是建模的基石。它告诉我们，杆柱上任意一点的加速度（方程左边）由两部分决定：一是该点附近杆段的应力差导致的恢复力（a² * ∂²u/∂x²），二是运动过程中受到的粘滞阻尼力（-c * ∂u/∂t）。理解每一项的物理意义，是后续设置正确边界条件和参数的关键。

2.2 有限差分法：将连续方程转化为MATLAB可解的代数方程

波动方程是偏微分方程（PDE），我们需要将其离散化才能用计算机求解。有限差分法（FDM）是最直观和常用的方法。其核心思想是用离散的网格点来逼近连续的空间和时间域。

1. 空间离散：将长度为L的杆柱等分为N段，空间步长Δx = L/N。离散点编号为i = 0, 1, 2, ..., N，其中i=0对应井口（光杆），i=N对应井下泵处。
2. 时间离散：将时间等分，时间步长为Δt。时间层编号为k = 0, 1, 2, ...。
3. 差分格式：我们采用中心差分格式来近似方程中的偏导数。这是精度和稳定性的一种平衡选择。
   • 时间二阶导数：∂²u/∂t² ≈ (u(i, k+1) - 2u(i, k) + u(i, k-1)) / (Δt)²
   • 空间二阶导数：∂²u/∂x² ≈ (u(i+1, k) - 2u(i, k) + u(i-1, k)) / (Δx)²
   • 时间一阶导数（阻尼项）：∂u/∂t ≈ (u(i, k+1) - u(i, k-1)) / (2Δt)

将上述差分格式代入波动方程，经过整理，我们可以得到关于未来时刻k+1层位移u(i, k+1)的显式递推公式：

u(i, k+1) = [ (a²Δt²/Δx²) * (u(i+1, k) + u(i-1, k)) + (2 - 2*a²Δt²/Δx² - cΔt) * u(i, k) + (cΔt/2 - 1) * u(i, k-1) ] / (1 + cΔt/2)

这个公式是MATLAB迭代计算的核心。它表明，杆上某一点下一时刻的位移，可以由当前时刻及前一时刻该点及其相邻两点的位移计算出来。

2.3 稳定性条件：CFL数的关键作用

显式差分格式是有条件稳定的。稳定性条件由著名的CFL（Courant-Friedrichs-Lewy）条件决定：

a * Δt / Δx ≤ 1

这个条件有深刻的物理意义：数值计算中信息传播的速度（Δx/Δt）必须大于或等于物理世界中波传播的实际速度（a）。否则，物理信息还未传到，计算就已经进行了，必然导致结果发散（数值爆炸）。在实际编程中，我们通常取CFL = a * Δt / Δx = 0.8 ~ 0.95，在保证稳定的前提下获得较大的时间步长以提高计算效率。

在MATLAB中，我们需要先根据杆长L、波速a和期望的空间分辨率（N）来确定Δx，然后根据CFL条件反推出最大允许的Δt。这是仿真能否成功运行的第一个关键检查点。我曾因为初期忽略了此条件，设置了过大的Δt，导致计算几十步后位移值就溢出成Inf或NaN，排查了许久才定位到这个基础问题。

3. 边界条件与载荷：连接模型与现实的桥梁

离散化的方程描述了杆柱内部的动力学，而系统的行为最终由边界（井口和泵端）驱动。边界条件的设置直接决定了模型的输入，是仿真是否“真实”的另一个生命线。

3.1 地面边界条件：驴头运动规律

在井口（i=0），位移u(0, t)是已知的，它由抽油机的几何结构和电机运动决定。最常见的简化是将其视为简谐运动：

u(0, t) = 0.5 * S * [1 - cos(2π * t / T)]

其中S是冲程（米），T是冲程周期（秒）。在MATLAB实现时，我们直接在每一个时间步k，将上述公式计算出的值赋给u(0, k)。更精确的模型可以考虑游梁式抽油机的实际连杆机构运动，用三角函数组合来描述，但简谐运动在多数诊断场景下已足够。

3.2 泵端边界条件：流体载荷与泵阀动力学

泵端（i=N）的边界条件最为复杂，因为它涉及到杆柱与井下流体、泵阀的相互作用。这是建模的难点和重点。我们通常采用力平衡条件。

在泵处，杆柱的力（由胡克定律计算：F_rod = EA * ∂u/∂x |_{x=L}）必须与作用在泵上的流体载荷F_fluid平衡。

F_rod = F_fluid

流体载荷F_fluid由泵筒内的压力决定：F_fluid = A_plunger * (P_below - P_above)

其中：

• A_plunger是柱塞截面积。
• P_below是泵吸入口压力（与地层流压相关）。
• P_above是泵排出口压力（与油管液柱压力、井口回压相关）。

而泵阀的开启与关闭，会动态改变P_above和P_below。一个相对实用且计算量可接受的简化模型是：

1. 上冲程：柱塞上行。固定阀打开，游动阀关闭。P_below≈ 吸入口压力（较低），P_above为上一冲程排入油管的液柱压力。此时F_fluid方向向上，是最大载荷。
2. 下冲程：柱塞下行。固定阀关闭，游动阀打开。P_above≈P_below≈ 吸入口压力。此时F_fluid很小甚至为负（向下），是最小载荷。

在MATLAB代码中，我们需要在每个时间步判断泵阀状态，从而动态计算F_fluid，然后利用力平衡条件推导出泵端 (i=N) 的位移u(N, k+1)。这通常需要将泵端的力平衡方程与内部节点的差分方程联立求解，或通过迭代逼近。

实操心得：泵阀逻辑的实现是代码中最易出错的部分。一个常见的错误是阀的开关判断逻辑与杆柱运动速度∂u/∂t的符号未正确关联，导致载荷曲线出现违背物理规律的震荡。我的调试方法是：单独输出一个冲程周期内泵端的位移、速度、计算出的上下压力及阀状态，绘制在一张图上，人工检查其逻辑时序是否正确。这个过程虽然繁琐，但一劳永逸。

4. MATLAB仿真实现：从零构建诊断原型

有了理论框架，我们开始在MATLAB中将其实现。我将以构建一个基础仿真原型为例，详解关键步骤和代码片段。

4.1 参数初始化与网格生成

% 1. 系统基本参数 L = 1000; % 杆柱长度，米 E = 2.1e11; % 钢弹性模量，Pa rho = 7850; % 钢密度，kg/m^3 a = sqrt(E/rho); % 波速，m/s c = 0.5; % 经验阻尼系数，1/s S = 3.0; % 冲程，m T = 10; % 周期，s f = 1/T; % 频率，Hz omega = 2*pi*f; % 角频率，rad/s % 2. 离散化参数 N = 100; % 空间分段数 dx = L / N; % 空间步长 CFL = 0.9; % CFL数，取0.9保证稳定 dt = CFL * dx / a; % 由稳定性条件确定时间步长 total_time = 5 * T; % 仿真总时间，模拟5个冲程 M = round(total_time / dt); % 总时间步数 % 3. 初始化位移场 (空间点数 N+1, 时间步数 M+1) u = zeros(N+1, M+1); % u(i, k), i:空间索引， k:时间索引

4.2 核心迭代循环与边界处理

% 4. 初始条件（假设从静止开始，杆柱处于拉伸平衡位置） % u(:, 1) 和 u(:, 2) 可设为相同值，或根据静载荷计算一个初始位移分布 % 这里简单设为0 u(:, 1) = 0; u(:, 2) = 0; % 预计算系数，提高循环效率 r = a * dt / dx; coeff1 = r^2; coeff2 = 2 - 2*r^2 - c*dt; coeff3 = c*dt/2 - 1; denom = 1 + c*dt/2; % 5. 时间迭代主循环 for k = 2:M % k代表当前已知的“过去”层，我们要计算 k+1 层 % 5.1 更新地面边界 (i=0) - 简谐运动 t = (k-1) * dt; % 当前时间 u(1, k+1) = 0.5 * S * (1 - cos(omega * t)); % 注意MATLAB索引从1开始 % 5.2 更新内部节点 (i=2 到 i=N) for i = 2:N u(i, k+1) = ( coeff1*(u(i+1, k) + u(i-1, k)) + ... coeff2*u(i, k) + ... coeff3*u(i, k-1) ) / denom; end % 5.3 更新泵端边界 (i=N+1) - 这里需要嵌入泵阀模型 % 这是一个简化示例，假设泵端力已知或通过简单关系给出 % 实际中这里需要调用一个独立的函数来计算泵端载荷和位移 [u(N+1, k+1), pump_force(k)] = updatePumpBoundary(u(N, k), u(N+1, k), u(N+1, k-1), dt, dx, E, A_rod, ...); end

updatePumpBoundary函数是工程实现的核心，它封装了第3.2节所述的泵阀逻辑和力平衡计算。由于其复杂性，它本身可能包含条件判断、压力计算和方程求解。

4.3 关键结果的可视化：诊断的“眼睛”

仿真完成后，我们必须将数据转化为直观的图形，这是诊断分析的基础。

% 6. 结果可视化 time_axis = (0:M) * dt; position_axis = (0:N) * dx; % 6.1 地面示功图（载荷-位移图） % 计算光杆载荷：F_surface = E * A_rod * (u(2,:) - u(1,:)) / dx （一阶差分近似应力） F_surface = E * A_rod * (u(2, :) - u(1, :)) / dx; figure; plot(u(1, :), F_surface / 1000, 'b-', 'LineWidth', 1.5); % 载荷单位化为kN xlabel('光杆位移 (m)'); ylabel('光杆载荷 (kN)'); title('仿真地面示功图'); grid on; % 6.2 泵端示功图 % 计算泵端载荷（已在循环中存储于 pump_force 数组） figure; plot(u(N+1, :), pump_force / 1000, 'r-', 'LineWidth', 1.5); xlabel('泵位移 (m)'); ylabel('泵载荷 (kN)'); title('仿真泵端示功图'); grid on; % 6.3 杆柱应力/位移分布动画 (可选，但非常直观) figure; for k = 1:50:M+1 % 每隔50帧显示一帧 plot(position_axis, u(:, k), 'b-o'); xlabel('杆柱位置 (m)'); ylabel('位移 (m)'); title(['杆柱位移分布 @ t = ', num2str((k-1)*dt, '%.2f'), ' s']); ylim([-0.5*S, 1.5*S]); grid on; drawnow; pause(0.05); end

将仿真得到的地面示功图与现场实测的示功图进行对比，是验证模型准确性的第一步，也是故障诊断的起点。

5. 基于模型仿真的故障诊断方法论

建立了准确的模型后，我们就可以将其用作一个“数字孪生体”。诊断的基本思路是：对比分析。

5.1 故障注入与特征库构建

健康的系统模型在标准参数下运行，会生成一个“基准”示功图。当某种故障发生时（如泵漏失、杆柱断脱、油管锚失效等），系统的某个或某几个参数会发生变化。我们在仿真模型中人为地修改这些参数，模拟故障状态，并记录下对应的示功图变化。

在MATLAB中，我们可以编写一个脚本，循环遍历多种故障模式及其严重程度，批量运行仿真，自动提取每个故障示功图的特征参数（如最大载荷、最小载荷、图形面积、斜率变化点等），构建一个“故障模式-特征向量”数据库。

5.2 实测数据对比与智能诊断

当获得一口井的实测地面示功图后，诊断流程如下：

1. 数据预处理：使用MATLAB的smoothdata函数滤除噪声，用findpeaks函数定位冲程起点和终点，进行归一化对齐。
2. 特征提取：从预处理后的实测图中提取与仿真特征库相同的一组特征参数。
3. 相似度匹配：将实测特征向量与故障特征库中的每一个向量进行相似度计算。简单的方法可以用欧氏距离，复杂一些可以用动态时间规整（DTW）来比较整个图形序列，MATLAB有dtw函数。
4. 故障识别与置信度评估：找到相似度最高的故障模式，并计算其匹配置信度。可以设定一个阈值，低于阈值则判定为“未知故障”或“复合故障”。

% 简化版诊断匹配示例 % measured_features: 从实测图提取的1xm特征向量 % fault_library: n x m 矩阵，n种故障，每种有m个特征 % fault_labels: n x 1 细胞数组，存储故障名称 distances = zeros(size(fault_library, 1), 1); for i = 1:size(fault_library, 1) distances(i) = norm(measured_features - fault_library(i, :)); end [min_dist, idx] = min(distances); confidence = 1 / (1 + min_dist); % 一个简单的置信度计算 if min_dist < threshold diagnosed_fault = fault_labels{idx}; fprintf('诊断结果：%s， 置信度：%.2f%%\n', diagnosed_fault, confidence*100); else fprintf('未匹配到已知故障模式。\n'); end

5.3 诊断系统的工程化考量

将上述原型发展为可用的诊断系统，还需要考虑以下工程问题：

• 模型参数标定：模型的准确性依赖于输入参数（如阻尼系数c、泵效、流体性质等）。我们需要利用井史数据或某次健康的实测数据，通过反演算法（如MATLAB的fminsearch,lsqnonlin）来校准这些参数，使仿真示功图与健康状态实测图最佳匹配。这是一个“模型调参”的过程。
• 实时性要求：对于在线监测，仿真速度必须快。可以考虑使用更高效的数值方法（如特征线法），或采用降阶模型（ROM）。MATLAB Coder可以将核心仿真代码转换为C/C++，显著提升速度。
• 不确定性处理：实测数据噪声大，模型本身也有简化误差。诊断结果应给出概率或置信区间，而不是绝对断言。可以引入贝叶斯推理框架来融合多源不确定信息。

在我参与的一个项目中，我们利用校准后的模型，成功预警了一口井的杆柱偏磨加剧趋势。模型仿真显示，在特定位置杆柱的横向振动加剧，对应应力集中。我们建议调整了扶正器位置，避免了后续可能发生的杆断事故。这种从“事后诊断”到“事前预警”的转变，正是数学建模价值的最高体现。

6. 超越基础：模型进阶与集成应用

基础的一维波动方程模型已经能解决大部分问题，但对于更复杂的场景，我们需要对其进行扩展。

6.1 三维振动与屈曲分析

在斜井、水平井或存在严重摩阻的情况下，杆柱不仅纵向振动，还会发生横向振动和螺旋屈曲。这就需要建立三维的杆柱动力学模型，控制方程将扩展到一组耦合的偏微分方程。虽然计算量剧增，但MATLAB的PDE Toolbox或自己编写有限元代码可以应对。这类模型可以精确分析杆管偏磨位置，为优化扶正器布置提供定量依据。

6.2 与地下渗流模型耦合

抽油系统的效率最终受制于地层供液能力。我们可以将井筒杆柱模型与描述地层流体向井筒流动的渗流模型（如使用达西定律）耦合起来。泵吸入口的压力P_below不再是一个固定值或简单假设，而是由地层渗透率、流体粘度、生产压差等参数动态计算得出。这种耦合仿真能更真实地模拟“供液不足”等生产动态问题。

6.3 集成到SCADA与数字孪生平台

最终的愿景是将这个MATLAB诊断模型集成到油田的SCADA（数据采集与监控系统）或数字孪生平台中。实现路径可以是：

1. MATLAB Production Server：将诊断算法部署为Web服务。SCADA系统实时将示功图数据发送到该服务，并接收返回的诊断结果和健康评分。
2. 编译独立应用：使用MATLAB Compiler将整个诊断程序打包成独立的可执行文件（.exe）或动态链接库（.dll），供其他平台（如C#、Java开发的平台）调用。
3. 云化部署：结合MATLAB Online或第三方云平台，实现多口井数据的集中化、规模化诊断分析。

从一行行数学公式和MATLAB代码，到一个能够精准模拟物理现实、并用于指导生产的诊断系统，这个过程充满了挑战，但也极具成就感。它要求我们不仅懂数值计算和编程，更要深入理解抽油系统背后的每一个物理细节。每一次模型的调优、每一次诊断的成功验证，都是对“机理+数据”这一工业智能化路径的坚实印证。当你看到自己构建的模型输出的示功图，与现场高精度传感器录制的图形高度重合时，那种跨越虚拟与现实的连接感，正是工程建模的魅力所在。