企业官网建设流程全解析

用MATLAB做三轴试验强度包线：最小二乘法拟合φ和c的完整套路

做三轴剪切试验写报告或者做毕设的同学，到了数据处理这一关，十有八九会被“画莫尔圆、作公切线、读内摩擦角φ和粘聚力c”这套流程折磨过。早几年我自己也是拿直尺和三角板在坐标纸上一条条贴公切线，贴完还要用两把尺子量角度，费时不说，不同人贴出来的切线位置差个一两度都是常事——而这一两度换算到φ和c上面，结论可能就差出一大截。后来我彻底改成用MATLAB处理：输入每级围压和破坏时主应力差，用最小二乘法自动拟合强度包线，φ和c直接出数值，图也顺手画好，省下来的时间用来复查数据质量，比手工画图划算得多。这篇就把整套思路和代码完整写出来，讲清楚每一行在干什么、为什么这么算，以及哪些坑我踩过之后希望你直接绕开。

先说适用对象：手里有三轴试验数据、需要按摩尔-库仑准则求总应力或有效应力抗剪强度指标（内摩擦角φ、粘聚力c）的本科生、研究生和工程师。代码只用到MATLAB最基础的polyfit和plot函数，不依赖任何额外工具箱，R2016b之后的版本都能跑。即使你MATLAB还不熟，按下面的步骤抄也能出图出数。

1. 先搞懂原理：三轴试验、莫尔圆和强度包线怎么串起来

1.1 三轴试验到底能给我们哪些数

常规三轴压缩试验的做法，是把3到4个相同土样分别放在不同的恒定围压σ3下，然后施加轴向压力让试样剪切破坏。围压就是压力室里的水压，一般取100、200、300、400 kPa这一档。每个试样剪切过程中，仪器记录轴向应变和主应力差（σ1-σ3）的关系曲线，试验结束时从曲线上取破坏点对应的主应力差（σ1-σ3）f——如果曲线有峰值就取峰值，没有峰值就按15%轴向应变对应的值取。

有了围压σ3和破坏时主应力差（σ1-σ3）f，就能算出试样破坏时的大主应力：

σ1f = σ3 + (σ1-σ3)f

一组围压对应一个试样，一个试样就能画出一个莫尔圆。比如我用过的典型数据：

这就是全部输入数据了。后面所有的圆、所有的拟合直线，都是从这四个数推出来的。

1.2 莫尔圆和强度包线究竟是什么关系

一个试样破坏时，破坏面上的正应力和剪应力并不是任意值，而是落在莫尔圆上。莫尔圆的圆心在σ轴上，横坐标是p = (σ1f + σ3) / 2，半径是q = (σ1f - σ3) / 2。在τ-σ坐标系里画这个圆，圆上每一个点就代表该试样在某一个方向平面上的应力状态；圆顶点的含义就是最大剪应力对应的面，也就是与大主应力面夹角45°的平面。

摩尔-库仑强度理论的核心假设是：土体破坏时，强度包线是τ = c + σ·tanφ这条直线，其中c是粘聚力，φ是内摩擦角。当某个莫尔圆恰好和这条包线相切时，土体就处于极限平衡状态。所以数据处理的目标很明确：找到一条直线，使之与所有试验得到的莫尔圆都“尽可能相切”，然后读出这条直线的截距和倾角，就是c和φ。

1.3 别再用直尺贴公切线了

手画公切线的经典流程是：把所有莫尔圆画在同一张坐标纸上，用直尺在圆族外侧找一条能与所有圆相切的直线，然后用量角器量倾角，读纵截距。这个方法最大的问题是主观性太强。实际试验数据很少能完美落在同一条切线上，总有一个圆稍微凸出一点，另一个圆稍微缩进去一点，直尺到底贴哪一个，不同人贴法不同，结果也就不同。尤其当某个圆的数据质量不佳时，人眼容易被异常点带偏。

最小二乘法做的事情，本质上就是用“残差平方和最小”这个明确标准，替代人眼的直觉判断。拟合结果不受操作者影响，同一份数据谁跑都一样，还能顺手给出拟合优度指标，用来判断整组试验数据靠不靠谱。这就是我推荐用MATLAB跑最小二乘而不是手工贴线的原因。

2. 拟合方法选型：把“找公切线”变成“拟合一条直线”

2.1 关键的坐标变换：从莫尔圆到q-p平面

直接对莫尔圆族求公切线在数学上是个非线性问题，新手很容易卡在这里。工程上更常用的做法是绕一步：不直接拟合包线，而是对每个莫尔圆的圆心横坐标p和半径q做线性回归，拟合这条线叫Kf线，在q-p坐标系里画出来就是一条直线。这条线贯穿各个莫尔圆的顶点，所以有些教材也叫它“顶点线”或“强度线”。

为什么可以绕这一步？回到相切条件。设强度包线为τ = c + σ·tanφ，任意一个莫尔圆圆心在(p, 0)，半径为q。圆心到包线的距离等于半径时，莫尔圆与包线相切。把点到直线的距离公式代进去整理，会得到一个极其干净的线性关系：

q = c·cosφ + p·sinφ

也就是说，在q-p平面上，各莫尔圆的半径q和圆心横坐标p之间本身就呈线性关系，斜率是sinφ，截距是c·cosφ。所以我们对(p, q)数据点做一元线性回归，再把斜率、截距反算回φ和c，就等价于找到了一条与所有圆“相切程度最优”的强度包线。

2.2 斜率、截距怎么换算成φ和c

设一元回归得到的直线是：

q = a + k·p

对照上面的公式：

k = sinφ，a = c·cosφ

于是：

φ = arcsin(k)

c = a / cosφ

举我前面那组数据为例。算出的4个莫尔圆圆心横坐标和半径分别是：

对这4个点做最小二乘线性拟合，斜率和截距大约是k≈0.662、a≈13.6。然后换算：

• φ = arcsin(0.662) ≈ 41.4°
• cos(41.4°) ≈ 0.75，所以c = 13.6 / 0.75 ≈ 18.1 kPa

这是一组典型的“有一定粘聚力、内摩擦角较大”的土体参数。如果斜率偏小，φ就小，说明土的摩擦强度低；如果截距偏大，c就大，说明土体本身具备较好的“黏性”强度。

这里有个细节值得强调：三角函数计算时，MATLAB默认角度单位都是弧度。asin返回的是弧度，必须乘以180/pi才是度数；而算c = a / cos(phi_rad)时必须用cos(弧度值)，不要先转成角度再算cos，那样会得到完全错误的结果。我第一次写这个脚本就在这个细节上翻过车。

2.3 什么时候需要考虑更高级的拟合方法

普通最小二乘对三轴试验数据基本够用，但有一个前提：各级围压下的数据误差相对均匀。如果你发现低围压的数据点明显更离散、高围压点更集中，可以考虑给每个点分配不同权重，做加权最小二乘；如果有个别异常点严重偏离趋势，也可以用稳健回归减少它的影响。不过这两种方法对常规课程设计和大多数科研项目来说属于“升级选项”，不是必选项。我自己的习惯是先用普通最小二乘出结果，然后看R²和残差，如果R²掉到0.95以下再回头检查是不是有数据点本身有问题，而不是急着换高档算法。

3. MATLAB完整实现：从试验数据到φ、c和包线图

3.1 数据准备与起步代码

先把输入数据准备好。最直接的方式是在代码开头用数组把围压和主应力差写死，数据量不多的时候这样最清晰。如果数据量很大，建议放到Excel里用readmatrix读入，但核心计算逻辑完全一样。

% 三轴试验原始数据，单位统一为 kPa sigma3 = [100; 200; 300; 400]; % 各级围压 dsigma = [285; 455; 625; 790]; % 破坏时主应力差 % 计算破坏大主应力、莫尔圆圆心横坐标p与半径q sigma1 = sigma3 + dsigma; p = (sigma1 + sigma3) / 2; % 圆心横坐标 q = (sigma1 - sigma3) / 2; % 莫尔圆半径

注意单位必须统一。如果你实验室仪器读数是MPa，就全部用MPa；是kPa就全部用kPa，中途尽量不要混着写，否则算出来的c和φ数值会很离谱。

3.2 最小二乘拟合与指标换算核心代码

最小二乘线性拟合，最省事的实现是polyfit。polyfit(x, y, 1)返回的是降幂排列的系数向量，第一个元素是斜率，第二个元素是截距：

coef = polyfit(p, q, 1); k = coef(1); % 斜率，即 sin(phi) a = coef(2); % 截距，即 c * cos(phi) % 换算抗剪强度指标 phi_rad = asin(k); % 内摩擦角，弧度 phi_deg = phi_rad * 180 / pi; % 内摩擦角，角度 c = a / cos(phi_rad); % 粘聚力 % 拟合优度 R^2 q_fit = polyval(coef, p); res = q - q_fit; SS_res = sum(res.^2); SS_tot = sum((q - mean(q)).^2); R2 = 1 - SS_res / SS_tot; fprintf('内摩擦角 phi = %.2f°\n', phi_deg); fprintf('粘聚力 c = %.2f kPa\n', c); fprintf('拟合优度 R^2 = %.4f\n', R2);

这里有几个点值得展开说一下：

1. polyfit本身不需要统计工具箱，属于MATLAB基础函数，所以哪怕你用的是学校机房精简版也没关系。
2. R²的计算方法就是统计学里那个标准定义：1减去残差平方和与总平方和的比值。R²越接近1，说明(p, q)数据点越接近一条直线，也就是各莫尔圆的尺寸和位置关系越符合摩尔-库仑直线包线的假设。
3. 如果数据量少，比如只有2个围压点，R²必然等于1，但这不代表结果可靠，只是两点决定一条直线的必然结果。至少3到4级围压的R²才有评价意义。

3.3 绘制莫尔圆、强度包线与切点标注

绘图是整个流程里比较能体现细节的部分。首先，绘制莫尔圆推荐用参数方程：一个圆心在(p, 0)、半径为q的圆，参数方程可以写成σ = p + q·cos(2θ)，τ = q·sin(2θ)，其中θ对应实际物理面的法线与大主应力方向的夹角。θ从0取到π，画的就是上半圆；想画完整圆把θ改成0到2π即可。教科书里强度包线讨论的通常都是上半平面，所以我习惯只画上半圆，图面也干净。

figure('Color','w','Position',[100 100 800 600]); hold on; axis equal; grid on; xlabel('正应力 \sigma (kPa)'); ylabel('剪应力 \tau (kPa)'); title(sprintf('三轴试验莫尔圆与强度包线 \\phi=%.1f°, c=%.1f kPa', phi_deg, c)); colors = lines(length(p)); % 为每个圆分配不同颜色 for i = 1:length(p) theta = linspace(0, pi, 200); sigma_circle = p(i) + q(i) * cos(2 * theta); tau_circle = q(i) * sin(2 * theta); plot(sigma_circle, tau_circle, 'Color', colors(i,:), 'LineWidth', 1.5); % 标记莫尔圆顶点 plot(p(i), q(i), 'o', 'MarkerFaceColor', colors(i,:), 'MarkerEdgeColor', 'k'); end % 绘制真实强度包线 tau = c + sigma * tan(phi) sigma_line = linspace(0, max(p) + max(q), 200); tau_line = c + sigma_line * tan(phi_rad); plot(sigma_line, tau_line, 'k-', 'LineWidth', 2.5);

这段代码有3个细节容易踩坑：

第一，axis equal必须加。不加的话，图形窗口会根据σ和τ的取值范围自动拉伸坐标轴，明明应该是圆形的莫尔圆会被画成椭圆。这个坑我见很多同学踩过，圆一变形，包线和圆是否相切就完全看不出来了。

第二，绘制包线用的是tan(phi_rad)，角度单位必须是弧度。如果你存了phi_deg就直接写tan(phi_deg)，算出来的斜率完全不对，画出来的包线会飘到天上去。

第三，不要在莫尔圆图上直接把Kf线画出来。Kf线是在q-p坐标系里的直线，坐标轴是p和q；σ-τ坐标系里画莫尔圆时，如果把Kf线原样叠上去，它的位置并不对应强度包线，很容易给读者或者审稿人造成误解。想展示拟合效果，可以单独画一张q-p散点加拟合直线的图；在主图上只画真实包线τ = c + σ·tanφ就够了。

另外，如果想标注强度包线与每个莫尔圆的切点，公式也不复杂。切点位于上半圆，对应中心角2θ = 90° + φ，所以切点坐标是：

σ* = p - q·sinφ τ* = q·cosφ

把这个点画在每个莫尔圆上，可以很直观地看到包线是否真的与圆相切：

% 计算并标记切点位置 for i = 1:length(p) sigma_t = p(i) - q(i) * sin(phi_rad); tau_t = q(i) * cos(phi_rad); plot(sigma_t, tau_t, 'rx', 'MarkerSize', 8, 'LineWidth', 1.5); end

如果拟合效果好，这些叉号会正好落在黑色包线上；如果某个点偏离明显，说明该级围压的试验数据可能有异常，值得回去查原始记录。

3.4 输出结果说明与拟合效果判断

运行完整脚本后，命令行会输出类似这样的结果：

内摩擦角 phi = 41.42° 粘聚力 c = 18.10 kPa 拟合优度 R^2 = 0.9993

R²达到0.9993说明这组数据线性程度很好，σ3取100到400 kPa范围内，摩尔-库仑直线包线的假设是合适的。如果R²偏低，要分情况看待：一是数据本身波动大，二是强度包线在这个应力范围内本身就明显弯曲。对于后者，比如某些超固结黏土在低围压段的包线呈现明显曲率，这时候强行用直线拟合虽然也能出数，但要把适用应力范围在报告里写清楚，不能外推到围压范围之外。

4. 实际操作中的常见坑与排查建议

4.1 莫尔圆被画成“椭圆”，八成是坐标轴比例问题

sigma轴和tau轴虽然单位都是kPa，但取值范围不同：sigma从0到1000多，tau最大也就400多。MATLAB默认会按数据范围自动缩放坐标轴，两个轴的单位长度不一致，圆立刻变扁或者变瘦。解决方式就是前面代码里那句axis equal。加了之后，两个轴按相同比例缩放，圆才是真正的圆。这个坑非常好避免，但几乎每次有人把图发给我看，第一眼就能看到圆变形——所以无论如何，画莫尔圆之前先把这句话写上。

4.2 拟合斜率超过1怎么办

理论上k = sinφ，打死也不可能大于1。但试验数据是有噪声的，当某个围压级别数据偏差太大时，最小二乘拟合出来的斜率完全可能变成1.02甚至更大。这时候asin(k)会返回复数，或者直接报错。遇到这种情况不要慌，先按以下顺序排查：

1. 检查数据录入是否有误，比如把395错录成935。
2. 检查该级围压试样的破坏模式。如果试样出现明显端部约束或者先沿着某个软弱面破坏，数据点会系统性偏离趋势线。
3. 检查是不是不同试样之间土体不均匀。三轴试验要求同一组试样初始状态尽量一致，如果孔隙比相差很大，强度自然对不上。

如果所有检查都找不到明确错误，可以作一个“手工约束版”的拟合，强制斜率不大于0.99，同时输出警告。但说实话，这种情况更可能的结论是这组试验不宜拿来求抗剪强度指标，重做试样可能比强行拟合更有意义。

4.3 总应力指标与有效应力指标千万别混用

UU、CU、CD三类三轴试验得到的c、φ含义完全不同。UU试验得到的是不排水总应力指标，适用于模拟饱和黏性土快速加载的短期稳定问题；CU试验如果不同步测孔压，得到的也是总应力指标，通常记为ccu和φcu；CD试验以及CU试验中同步测定孔压后换算出的结果是有效应力指标，记为c'和φ'，适用于排水长期工况。画强度包线时，同一张图里要么全是总应力莫尔圆，要么全是有效应力莫尔圆，绝对不能混用。我见过一份课程报告，低围压用有效应力圆、高围压用总应力圆，拟合出来的φ和c完全不知道属于什么工况，这种数据拿到答辩现场基本就是送分题给老师挑毛病。

4.4 数据点偏少、围压范围偏窄时怎么处理

如果只有两个围压点，比如σ3 = 100和300 kPa，那么两个点一定能连出一条直线，φ和c看起来也有模有样，但这只是被两个点“硬锁”出来的结果，没有任何冗余来验证合理性。常规三轴试验都要求至少3到4级围压，并且围压范围要覆盖你关心的应力区间。整理报告时如果条件确实只允许做2个有效试样，建议在结论部分明确写明“本次仅由2级围压数据拟合，强度指标供参考”，避免给后续使用造成误导。

4.5 数据录入与单位引发的怪异结果

有一种很隐蔽的错误是把某组主应力差输成主应力。比如把σ1f = 385输入到dsigma的位置，相当于把莫尔圆半径翻了一倍，拟合出来的φ会异常偏高，c甚至会算出负值。碰到φ超过50°或者c出现负数的“可疑结果”，第一时间核原始数据，不要急着怀疑算法。我在处理同学的数据时，至少有一半的“异常结果”最后都出在录入错误上。

5. 批量处理和扩展方向

5.1 把拟合过程封装成函数，多组土样一次跑完

如果你要对多组土样分别求c和φ，最清晰的组织方式是把核心计算封装成一个函数，输入sigma3和dsigma，输出phi_deg、c和R2。这样主脚本只需要循环调用函数，读不同的数据文件，最后汇总结果到一个表格。

function [phi_deg, c, R2] = triaxial_fit(sigma3, dsigma) sigma1 = sigma3 + dsigma; p = (sigma1 + sigma3) / 2; q = (sigma1 - sigma3) / 2; coef = polyfit(p, q, 1); k = coef(1); a = coef(2); phi_rad = asin(k); phi_deg = phi_rad * 180 / pi; c = a / cos(phi_rad); q_fit = polyval(coef, p); SS_res = sum((q - q_fit).^2); SS_tot = sum((q - mean(q)).^2); R2 = 1 - SS_res / SS_tot; end

调用的时候，先读Excel，再循环：

data = readmatrix('试验数据.xlsx'); % 每一列依次是围压和主应力差 for g = 1:size(data, 3) % 按实际情况调整维度 [phi(g), c(g), R2(g)] = triaxial_fit(data(:,1,g), data(:,2,g)); end

这样几十组土样也能几分钟内出齐表格，比手动画图快一个量级。

5.2 从总应力指标扩展到有效应力指标

很多三轴试验课程项目做到CU试验时会记录孔隙水压力u。这时候可以分别用总应力(σ1, σ3)和有效应力(σ1' = σ1 - u, σ3' = σ3 - u)各算一组c、φ。两者的物理含义完全不同，在报告里都值得给出来。计算有效应力指标时，只要替换sigma3和dsigma的输入逻辑，把孔压扣掉就行，核心拟合代码完全不用改。

5.3 非线性强度包线的处理思路

如果你的土样在较宽的围压范围内数据显示包线明显弯曲，单纯直线拟合就不够了。一个常用的折中做法是把围压范围分段，每段分别做直线拟合，分别读取c、φ；另一个做法是改用非线性准则描述，这就超出本文最小二乘直线拟合的范畴了。对于大多数课程设计和常规岩土工程问题，摩尔-库仑直线包线仍然是通行做法，重点是把拟合的适用应力范围交代清楚。

最后分享一个我自己的小习惯：每次拟合完，我都会顺手把每个莫尔圆圆心到包线的垂直距离打印出来，和半径放在一起对比。垂直距离由公式d = (p·tanφ + c) / sqrt(tan²φ + 1)算出，理论上应该等于半径q。如果两者差值都在0.5 kPa以内，说明这组拟合线的“相切程度”是可信的；一旦某个圆明显偏离，一定是那级试验数据本身有情况。这个检查动作花不了两分钟，但能让你的结论在答辩或报告审查时经得起追问，我自己已经把它列进每次数据处理的标准流程了。