多波束数据体积变化计算

这段代码的作用是读取并处理多波束测量数据，对两次测量数据进行匹配，计算挖沙区域的体积变化。具体包括：

读取两次测量数据。
对每次测量数据按照x和y坐标进行排序。
进行间隔采样。
循环匹配两次测量数据中相近的点，并计算它们之间的距离。
对匹配结果进行滤波，去掉z变化小于阈值的点。
计算挖沙区域的体积变化，即计算两次测量数据中相同位置处的差值，并乘以相邻采样点间的体积。
可视化匹配结果、z变化和体积分布。

代码首先使用'read_excel'函数读取两次测量数据，分别存储在'samp1'和'samp2'变量中。然后，使用'sort_xy'函数对每次测量数据按照x和y坐标进行排序，方便后续匹配。接着，代码进行间隔采样，减少计算量。

接下来，代码进入循环匹配部分。对于'samp1'中的每一个点，代码遍历'samp2'中的所有点，寻找x和y坐标相近的点。找到匹配的点后，代码计算它们之间的距离，并存储在'R1'矩阵中。

在匹配过程中，代码使用了一个阈值'x_g'和'y_g'，用来判断两个点是否相近。如果两个点的x坐标差小于'x_g'，且y坐标差小于'y_g'，则认为这两个点是匹配的。

匹配结束后，代码使用'smoothdata'函数对z坐标进行平滑滤波，减少噪声的影响。然后，代码计算两次测量数据中相同位置处的z坐标差值，并将其存储在'delt_z'变量中。

为了进一步去除噪声，代码使用了一个阈值'delt_z_g'，将z坐标差值小于'delt_z_g'的点去除。

最后，代码计算挖沙区域的体积变化。代码计算了匹配点之间的体积，并对所有体积进行求和，得到总的体积变化。同时，代码还可视化了匹配结果、z变化和体积分布。

这段代码使用了多种数据处理和可视化方法，可以有效地计算挖沙区域的体积变化，并对结果进行分析和可视化。

clear
close all
addpath '多波束数据'

samp1=read_excel(2,1,1,'第一次 间距5m','全测');
samp2=read_excel(2,1,1,'第一次 间距5m','挖沙区域');

save samp.mat
load samp.mat
disp('读完文件');

  %对x y排序
  disp('对xy排序');
 [samp1] = sort_xy(samp1);
 [samp2] = sort_xy(samp2);
   x1=samp1(:,1);  y1=samp1(:,2);   z1=samp1(:,3);
   x2=samp2(:,1);  y2=samp2(:,2);   z2=samp2(:,3);
  plot3_view(x1,y1,z1,  '', '', '' )  
  title('sheet1原图');
  plot3_view(x2,y2,z2,  '', '', '' )
  title('sheet2原图');
   
x_g=10; %如果x相差x_g就认为相等
y_g=10; %如果x相差y_g就认为相等
jiange=1; %xy的取样跨度
[m1,n1]=size(samp1); % x y z
[m2,n2]=size(samp2); % x y z
%间隔采样
p1=[1:jiange:m1];
p2=[1:jiange:m2];
samp1=samp1(p1,:);
samp2=samp2(p2,:);
[m1,n1]=size(samp1); % x y z
[m2,n2]=size(samp2); % x y z

disp('开始处理匹配点');
%循环匹配，找xy相近的
R1=zeros(m1,9); %存储结果的内存， x1 y1 x2 y2 z1 z2  第1次采样的samp中的序号  第2次采样的序号   距离
for i=1:m1
    p2=0;
    R2=zeros(m2,9); %存储结果的内存， x1 y1 x2 y2 z1 z2  第1次采样的samp中的序号  第2次采样的序号   距离
    R2(:,end)= R2(:,end)+inf;
    x1=samp1(i,1);  y1=samp1(i,2);   z1=samp1(i,3);
    for j=1:m2
        x2=samp2(j,1);  y2=samp2(j,2);   z2=samp2(j,3);
        if abs(x1-x2)< x_g & abs(y1-y2) < y_g
            p2=p2+1;
            d=sqrt( (x1-x2)^2+(y1-y2)^2 );
            R2(p2,:)=[x1  y1  x2 y2  z1 z2  i j d];
        end
    end
    if p2 > 0  % 未挖沙前的采样点，找到了挖沙后的相近的点
        [~,idx]=min(R2(:,9));
        R1(i,:)=R2(idx,:); %保存最小的
         if mod(i,10000)==0
         fprintf('原来的测量第%d万个测量点已找到%d个相近的点
',i/10000,p2);
        end
    end
    
end

idx=find( R1(:,1)~=0);
R=R1(idx,:);

z1=R(:,5);
z2=R(:,6);
%对z,平滑滤波
z1=smoothdata(z1);
z2=smoothdata(z2);
delt_z=abs(z1-z2); %z的变化

delt_z_g=1e-2; %如果z的变化 小于delt_z_g 就认为是噪声就去掉
idx=find(delt_z >delt_z_g);
R=R(idx,:);
delt_z=delt_z(idx);

%计算体积
x=R(:,1);
y=R(:,2);
delt_x=abs(x(2:end)-x(1:end-1));
delt_y=abs(y(2:end)-y(1:end-1));
delt_z=delt_z(2:end);
v=delt_x.*delt_y.*delt_z;
sum_v=sum(v);

disp('总体积')
sum_v
 
x1=R(:,1); y1=R(:,2); z1=R(:,5);
x2=R(:,3); y2=R(:,4); z2=R(:,6);
delt_z=abs(z1-z2);

plot3_view(x1,y1,z2,  '', '', '' );
  title('相同点的图');
plot3_view(x1,y1,delt_z,  '', '', '' );
  title('delta-z图');
  
  plot3_view(x1,y1,v,  '', '', '' );
  title('计算的体积分布图');
  
  v1=v;
 save xyzv1.mat x1 y1 v1