matlab最小二乘法，绘制强度包线，并求解土体三轴试验抗剪强度指标 内摩擦角，粘聚力 前沟通好，代码为 莫尔圆，强度包络线

算三轴剪莫尔圆手切线的噩梦是不是还留着大二大三实验室报告册的红叉——要么三个圆勉强切一条歪歪扭扭的，要么老师让“再微调下内摩擦角试试”“粘聚力别勾到莫尔圆里头去啊同学”？

没关系，今天掏出MATLAB这台“土工绘图计算器”，一键生成准确到小数点后两位的强度包线，顺便搞明白最小二乘法在这儿是怎么给我们当“客观线评人”的。

先把咱们实验室攒的（或者随便编一组像真的的）常规固结不排水CU试验数据摆出来吧：

% 围压sigma3，单位kPa，假设做了3组：100, 200, 400
sigma3 = [100, 200, 400];
% 对应的破坏主应力差sigma1-sigma3，单位kPa
sigma1_diff = [220, 380, 710];

哦对了，先回忆下土力学里莫尔圆和强度包线的基本逻辑：破坏时的莫尔圆圆心是（(sigma1+sigma3)/2, 0），半径是（sigma1-sigma3)/2，而库仑强度方程是tau = c + sigma*tanφ。我们要找的就是让所有破坏时的剪应力点和包线偏差最小的c和φ——这偏差的衡量就是最小二乘法啦。

接下来分三步走：画莫尔圆、最小二乘拟合包线、最后把圆和线捏在一起再标出来关键指标。

第一步，画三个干干净净的圆。MATLAB画圆不像Python还要生成一堆theta和坐标点，自带个viscircles函数，输入圆心半径就行，还能直接调线宽颜色，太懂我们这种赶报告的懒人了。

先把sigma1算出来吧，差sigma1_diff加上sigma3就是了：

sigma1 = sigma3 + sigma1_diff;
centers = [(sigma1 + sigma3)/2, zeros(size(sigma3))]; % 圆心列矩阵？哦viscircles要行还是列？试下行更方便
radii = (sigma1 - sigma3)/2;

先别急着画包线，试试先画圆撑撑场面，顺便把x轴y轴调得像样点——别让圆只占个小角落，也别把剪应力搞成负数：

figure('Color','w','Position',[100 100 800 600]); % 白色背景，窗口大一点看得清
viscircles(centers, radii, 'Color','#4a79a4','LineWidth',2,'LineStyle','-'); % 用个柔和的土蓝色，不是土黄色哦莫尔圆用蓝色经典
hold on;
axis equal; % 必须加！不然圆会变成椭圆！
xlim([0, max(sigma1)+50]); % x轴留一点余量放包线的延伸
ylim([0, max(radii)+50]); % y轴同理
xlabel('\sigma (kPa)','FontSize',12,'FontWeight','bold');
ylabel('\tau (kPa)','FontSize',12,'FontWeight','bold');
title('常规CU试验莫尔圆与最小二乘强度包线','FontSize',14,'FontWeight','bold');
grid on; grid minor; % 加个小网格方便看刻度

现在运行到这儿，你应该能看到三个排列得整整齐齐的椭圆……哦不，加了axis equal之后的三个蓝色莫尔圆啦。

第二步，核心来了——最小二乘法拟合库仑强度包线。

很多同学可能会直接用MATLAB自带的polyfit对剪应力tau和正应力sigma做一次拟合？但这里有个小坑：库仑方程是莫尔圆的切线方程，不是随便拿破坏面上的（sigmaf, tauf）点去拟直线！哦不对等下——哦也对？因为每个破坏莫尔圆和库仑包线的切点就是（sigmaf, tauf），那如果我们能写出切点坐标和c、φ的关系，是不是就可以构造误差函数用最小二乘解了？

先推导一下破坏面的正应力和剪应力公式（土力学课必背但现在不用背，直接写）：

sigma_f = (sigma1+sigma3)/2 + (sigma1-sigma3)/2 * sinφ

tau_f = (sigma1-sigma3)/2 * cosφ

同时库仑方程要求tauf = c + sigmaf * tanφ

那我们把前两个代入第三个，消掉sigmaf和tauf，整理一下（这个推导过程别写太长太论文，大概说一句“把切点代入库仑方程消掉坐标，得到每个莫尔圆对应的误差残差”）：

残差ei = [(sigma1i - sigma3i)/2 cosφ] - [c + ((sigma1i + sigma3i)/2 + (sigma1i - sigma3_i)/2 sinφ) * tanφ]

哎再把tanφ换成sinφ/cosφ，化简一下残差会变得巨简单！分母没了，全是线性的？不对等下是关于c和sinφ的线性项！哦对哦，库仑方程可以变形为（tauf cosφ - sigmaf sinφ）= c，而每个莫尔圆的“理论左边”刚好就是这个值对吧？那残差ei = c + sigmafi sinφ - taufi cosφ，但又因为sigmafi和taufi可以用圆心Oi=(Oix, 0)、半径Ri表示成Oix + Ri sinφ和R_i cosφ，代入残差的话：

ei = c + (Oix + Ri sinφ) sinφ - R_i cosφ cosφ

= c + Oix sinφ + R_i (sin²φ - cos²φ)

哦cos²减sin²是cos2φ，但这里好像sin²减cos²也没关系，不过我们要的是最小二乘，就是要让Σei²最小对吧？那把ei²展开，对c和sinφ求偏导令其等于0，就能得到线性方程组啦！

哦对哦！这里不用迭代法求非线性方程，因为我们可以设两个未知量：k1 = c，k2 = sinφ，然后把残差改写为关于k1和k2的线性函数！

等下再仔细推一遍线性化（这一步一定要讲清楚，不然代码就像天书）：

刚才的化简过的切点代入残差（哦应该反过来，库仑方程的左边和莫尔圆的“理论左边切点值”的差，不过平方的话正负无所谓），刚才得到的：

matlab最小二乘法，绘制强度包线，并求解土体三轴试验抗剪强度指标 内摩擦角，粘聚力 前沟通好，代码为 莫尔圆，强度包络线

c + Oix sinφ = Ri cos²φ - Ri sin²φ + Oix sinφ？不对不对刚才的代入好像出错了，换个更简单的线性化思路吧——土力学里常用的“莫尔圆切线最小二乘等价于对（sigma, tau）做加权线性回归？不，或者用圆心和半径的关系构造线性方程组。

另一个更直观的办法：每个莫尔圆到强度包线的垂直距离（沿半径方向，因为半径是垂直于切线的）的平方和最小！对，这个才是几何意义上的最小二乘拟合，比随便找点更准确，因为有时候我们不知道破坏面的具体位置，但知道圆和切线的几何关系。

库仑强度包线写成标准直线形式：Ax + By + C = 0，不过我们用剪应力tau在y轴，正应力sigma在x轴的话，就是：tanφ x - 1 y + c = 0，也就是A=tanφ，B=-1，C=c。

一个点(x0,y0)到这条直线的距离是|A x0 + B y0 + C| / sqrt(A²+B²)，那莫尔圆的圆心(Oix, 0)到包线的距离，就必须等于莫尔圆的半径R_i（因为是切线嘛）！

所以几何意义上的残差平方和就是Σ[ ( |A Oix + B*0 + C| / sqrt(A²+B²) ) - R_i ]²，不过因为A、B、C有固定比例（B=-1，sqrt(A²+B²)=sqrt(tan²φ+1)=1/cosφ，绝对值也可以去掉，因为围压是正的，包线在上面，所以分子肯定是正的），所以可以简化这个残差：

把A=tanφ=sinφ/cosφ，B=-1，C=c代入，分子是（sinφ/cosφ）Oix + c，分母是1/cosφ，所以分子除以分母就是sinφ Oix + c * cosφ，这个值必须等于R_i！

哦！我的天刚才绕了一大圈，现在这个式子太清爽了：每个莫尔圆对应的sinφOix + cosφc = R_i！

这是什么？这是关于未知量（c, sinφ, cosφ）的方程？不对哦有三个未知量但只有sin²φ + cos²φ=1的约束？没关系没关系，我们可以把它看作关于（c*cosφ, sinφ）的线性方程！

令：

u = c * cosφ

v = sinφ

那每个方程就变成：Oix v + u = Ri → u + Oix v = Ri

哇！完全是线性方程组！三个方程（如果做了三组试验）两个未知量u和v，刚好可以用最小二乘法求最优解！

这个线性化太妙了对吧？完全不用迭代，直接A\b或者pinv(A)*b就能搞定，MATLAB最喜欢干这个了。

那现在把这个线性方程组的系数矩阵A和常数项b写出来：

% 构造线性方程组的系数矩阵A：第一列全是1（对应u），第二列是圆心的x坐标O_i_x（对应v）
A = [ones(length(sigma3), 1), (sigma1 + sigma3)/2];
% 常数项b就是每个莫尔圆的半径R_i
b = (sigma1 - sigma3)/2;
% 最小二乘求解u和v：A*X = b，X = [u; v]
X = A \ b;
u = X(1);
v = X(2);

现在u和v有了，接下来反算我们要的c和φ：

v是sinφ，所以φ = arcsin(v)，记得转换成角度哦，土力学报告里没人用弧度；

u是c*cosφ，cosφ可以用sqrt(1 - v²)算（因为v是sinφ，而且内摩擦角肯定在0到90度之间，cosφ是正的），所以c = u / sqrt(1 - v²)。

% 反算内摩擦角φ（弧度转角度）
phi_rad = asin(v);
phi_deg = phi_rad * 180 / pi;
% 反算粘聚力c
cos_phi = sqrt(1 - v^2);
c = u / cos_phi;
% 打印一下结果，保留两位小数
fprintf('常规CU试验最小二乘拟合结果：\n');
fprintf('内摩擦角φ = %.2f°\n', phi_deg);
fprintf('粘聚力c = %.2f kPa\n', c);

现在运行到这儿，命令窗口应该会跳出咱们的指标啦，比如刚才编的那组数据，算出来大概φ=25度左右，c=15kPa左右？自己可以试试。

第三步，画强度包线，顺便标上指标。

强度包线是直线，我们可以取x轴的两个点来画，比如x1=0（这时候y1=c，刚好是粘聚力的截距！这个点一定要标出来！），x2=max(sigma1)+50（刚才x轴留的余量），对应的y2=c + x2*tan(phi_rad)。

% 取强度包线的两个端点
x_line = [0, max(sigma1)+50];
y_line = c + x_line * tan(phi_rad);
% 画强度包线，用红色虚线？或者红色实线加粗？红色经典，加粗更显眼
plot(x_line, y_line, 'Color','#d62728','LineWidth',2.5,'LineStyle','-');
% 标上粘聚力的截距点(0,c)，画个小圆圈
plot(0, c, 'o','Color','#d62728','MarkerSize',8,'MarkerFaceColor','#d62728');
text(0+20, c+10, sprintf('c=%.2f kPa',c),'FontSize',11,'FontWeight','bold','Color','#d62728');
% 再标上内摩擦角，画个小角度标记
% 先找x轴附近的一个点计算角度标记的位置
theta = phi_rad;
% 取标记的半径，比如max(radii)/4
r_theta = max(radii)/4;
% 角度标记的弧线
x_theta = 0:0.01:r_theta*cos(theta);
y_theta = x_theta * tan(theta);
% 画弧线
plot(x_theta, y_theta, 'Color','#d62728','LineWidth',1.5);
% 标角度的文字
text(r_theta*cos(theta/2)+10, r_theta*sin(theta/2)+10, sprintf('φ=%.2f°',phi_deg),'FontSize',11,'FontWeight','bold','Color','#d62728');
hold off;

现在把所有代码连起来运行，一张完美的、有三个蓝色莫尔圆、一条红色加粗强度包线、标了c和φ的土力学三轴试验图就出来啦！命令窗口也有准确的数值，再也不用怕手切线手滑了。

对了刚才那组编的数据连起来的完整代码可以直接复制用哦，当然可以换成你自己的CU、UU或者CD试验数据，逻辑都是一样的——只要有围压sigma3和破坏主应力sigma1（或者sigma1-sigma3）就行。

希望这篇小分享能帮到正在赶土工试验报告的同学，或者复习土力学的考研党～