实际地铁隧道断面水平直径编程实例

说明

截取环中部位置断面计算环水平直径，截取隧道中心线垂直的断面点云时，可 以直接按y值截取，统计单环最大、最小y值，取平均即为环中间位置；截取断面数据时按一定厚度（如3-5cm）截取，取出满足要求的点云后，去掉y 坐标，只保留x，z，相当于将去除的点云投影到平面上； 在平面上拟合圆，得到圆心坐标及半径，需要剔除非隧道壁上的点云，可用选权迭代法剔除，

以拟合圆心的z坐标为参照，在断面数据中搜索用于拟合的扫描点，可截取z坐 标的上下10cm共20cm的断面点数据，利用截取的扫描点数据，选取合适的函数（三次多项式），编程拟合曲 线；用拟合圆心z坐标的高度，代入拟合曲线，求取拟合位置x值；同理，求取隧道另一侧拟合位置；用上面求得的隧道两侧的拟合点，计算拟合水平直径值。

---

关键：

1：拟合圆：

function [xc,yc,R] = circfit(x,y)
    %圆拟合函数
    % x^2+y^2+a(1)*x+a(2)*y+a(3)=0
    n=length(x); xx=x.*x; yy=y.*y; xy=x.*y;
    A=[sum(x) sum(y) n;sum(xy) sum(yy)...
    sum(y);sum(xx) sum(xy) sum(x)];
    B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
    a=A\B;            %x = A\B 用来求解线性方程 A*x = B.  A 和 B 的行数一致.
    xc = -.5*a(1);
    yc = -.5*a(2);
    R = sqrt((a(1)^2+a(2)^2)/4-a(3));  
end 

 2：剔点：先拟合圆，然后计算误差，每次剔除误差超限的点

 3：数据集划分：训练集和测试集：randperm函数打乱

4：多项式拟合预测：polyfit、polyval

---

代码

clear all;clc;
%% 数据读取
DATA=load("169_pc.pts");
%% 断面截取
Y_MAX=max(DATA(:,2));
Y_MIN=min(DATA(:,2));
Y_mean=(Y_MIN+Y_MAX)/2;
Y_d=Y_mean-0.03;Y_u=Y_mean+0.03;
sum0=1;
for i=1:size(DATA,1)
    if (DATA(i,2)<Y_u)&&(DATA(i,2)>Y_d)
        X(sum0,1)=DATA(i,1);Z(sum0,1)=DATA(i,3);
        sum0=sum0+1;
    end
end
figure(1)
plot(X,Z,"b.");
xlabel("X (m)")
ylabel("Z (m)")
title('断面截取结果')
%% 剔点
%先剔除隧道底部点云
sum0=1;
for i=1:size(X,1)
    if Z(i,1)>0
        X1(sum0,1)=X(i,1);Z1(sum0,1)=Z(i,1);
        sum0=sum0+1;
    end
end
figure(2)
plot(X1,Z1,"b.");
xlabel("X (m)")
ylabel("Z (m)")
title("底部剔除后结果图")
%剔除隧道墙壁上的点
X_TRUE=X1;Z_TRUE=Z1;
error=5;j=1;
sum1=sum0-1;
while max(error)>0.5
    error=[];
    [X0,Z0,R] = circfit(X_TRUE(1:(sum0-1),j),Z_TRUE(1:(sum0-1),j));
    j=j+1;sum0=1;
    for i=1:sum1
        error(i,1)=abs(sqrt((X_TRUE(i,j-1)-X0)^2+(Z_TRUE(i,j-1)-Z0)^2)-R);
        if error(i,1)<0.08
            X_TRUE(sum0,j)=X_TRUE(i,j-1);
            Z_TRUE(sum0,j)=Z_TRUE(i,j-1);
            sum0=sum0+1;
        end
    end
    sum1=sum0-1;
end
figure(3)
plot(X_TRUE(1:(sum0-1),j),Z_TRUE(1:(sum0-1),j),"b.");
hold on
plot(X0,Z0,"r*")
aplha=0:pi/40:2*pi;
x=X0+R*cos(aplha);
z=Z0+R*sin(aplha);
plot(x,z,'-');
xlabel("X (m)")
ylabel("Z (m)")
title("剔除后结果图")
%% 截取z坐标的上下10cm共20cm的断面点数据
Z_u=Z0+0.1;Z_d=Z0-0.1;
sum2=1;
for i=1:size(X_TRUE(1:(sum0-1),j),1)
    if (Z_TRUE(i,j)<Z_u)&&(Z_TRUE(i,j)>Z_d)
        X_MID(sum2,1)=X_TRUE(i,j);Z_MID(sum2,1)=Z_TRUE(i,j);
        sum2=sum2+1;
    end
end
figure(4)
plot(X_MID,Z_MID,"r*");
hold on
plot(X_TRUE(1:(sum0-1),j),Z_TRUE(1:(sum0-1),j),"b.");
plot(X0,Z0,"ro")
plot(x,z,'-');
xlabel("X (m)")
ylabel("Z (m)")
title("截取圆心高度结果图")
%% 三次多项式拟合
num_l=1;num_r=1;
for i=1:size(X_MID,1)
    if X_MID(i,1)<X0
        X_L(num_l,1)=X_MID(i,1);Z_L(num_l,1)=Z_MID(i,1);
        num_l=num_l+1;
    else
        X_R(num_r,1)=X_MID(i,1);Z_R(num_r,1)=Z_MID(i,1);
        num_r=num_r+1;
    end
end
%划分训练集测试集
%随机划分
X_L= X_L(randperm(size(X_L,1)),:);
Z_L= Z_L(randperm(size(Z_L,1)),:);
X_R= X_R(randperm(size(X_R,1)),:);
Z_R= Z_R(randperm(size(Z_R,1)),:);
%按比例分
ind_L = round(0.7 *size(X_L,1));
ind_R = round(0.7 *size(X_R,1));
%训练集
TRAIN_X_L=X_L(1:ind_L,1);
TRAIN_Z_L=Z_L(1:ind_L,1);
TRAIN_X_R=X_R(1:ind_R,1);
TRAIN_Z_R=Z_R(1:ind_R,1);
%测试集
TEST_X_L=X_L(ind_L+1:end,1);
TEST_Z_L=Z_L(ind_L+1:end,1);
TEST_X_R=X_R(ind_R+1:end,1);
TEST_Z_R=Z_R(ind_R+1:end,1);
%左侧
p_L = polyfit(TRAIN_Z_L,TRAIN_X_L,2);
%内精度
X_L_1 = polyval(p_L,TRAIN_Z_L);
vx_L_1=X_L_1-TRAIN_X_L;
dx_L_1=sqrt(sum(vx_L_1.^2)/(size(vx_L_1,1)-6));
%外精度，训练集上误差
X_L_2 = polyval(p_L,TEST_Z_L);
vx_L_2=X_L_2-TEST_X_L;
dx_L_2=sqrt(sum(vx_L_2.^2)/(size(vx_L_2,1)-6));
%预测
R_L=polyval(p_L,Z0);
%右侧
p_R = polyfit(TRAIN_Z_R,TRAIN_X_R,2);
%内精度
X_R_1 = polyval(p_R,TRAIN_Z_R);
vx_R_1=X_R_1-TRAIN_X_R;
dx_R_1=sqrt(sum(vx_R_1.^2)/(size(vx_R_1,1)-6));
%外精度，训练集上误差
X_R_2 = polyval(p_R,TEST_Z_R);
vx_R_2=X_R_2-TEST_X_R;
dx_R_2=sqrt(sum(vx_R_2.^2)/(size(vx_R_2,1)-6));
%预测
R_R=polyval(p_R,Z0);
%直径
RR=R_R-R_L;
DR=RR-2*R;
%% 成图
figure(4)
plot(X,Z,"b.");
hold on
plot(X0,Z0,"r*")
plot(x,z,'-');
XX=[X0-R,X0+R];ZZ=[Z0,Z0];
plot(XX,ZZ,'-');
text(X0-2,Z0+0.5,'环号：169  水平直径：5.4151m')
xlabel("X (m)")
ylabel("Z (m)")
title("地铁隧道断面水平直径计算结果")