代码分享：面波数据快速成图

前言

目前，物探数据主要用surfer软件成图，surfer软件具有强大的插值和绘图功能，成图比较美观。但是，生产过程中大量的物探数据，依靠excel和surfer来成图耗费人力时间成本。本博文在MATLAB平台上开发了一套用于面波数据快速成图的小程序，仅供大家借鉴。

1、成图效果展示

1.1 原始图像

对面波数据采用geogiga软件反演，导出视横波数据，在matlab中编辑克里金插值算法的代码，输出图像。

1.2 高程转换

将地表GPS测量的高程，与面波探测的深度进行转换，得到真实的高程。

1.3 里程换算

将地表GPS测量的里程与高程，与面波探测的深度和水平距离进行换算，由于面波测点在地表不等距分布，因此里程也是不等间距分布，换算之后得到真实的高程与里程。

1.4 图像加工

为了得到比较美观的图像，在MATLAB中对图像进行加工。

2、数据读取与图像保存

2.1 读取面波视横波速度数据

选择数据文件夹。

% 读取面波数据
[FileName,PathName] = uigetfile('*.txt','请选择视横波速度文件',...
    'MultiSelect','on');
filename = strcat(PathName,FileName);
data = importdata(filename);
fprintf('
读取视横波速度完成!
请按任意键继续...
');

提取数据，自编函数dealMBdata。

% 初始参数设置
% 最大深度
depth_max = 80;
% 插值点数
num_points = 40;

% 面波数据预处理
[points,vs_value,xlocation] = dealMBdata(data);

2.2 数据与图像保存

% 保存数据
clear xx yy zz
xx = X1(:);
yy = Y1_new(:);
zz = YX(:);
C = [xx,yy,zz];
dlmwrite(strcat(PathName,'mianbo.dat'),C);
clear yy
yy = Y1_new(1,:);
high = [xa',yy'];
dlmwrite(strcat(PathName,'gaocheng.dat'),high);

3、自编函数

3.1 dealMBdata函数

function [points,vs_value,xlocation] = dealMBdata(data)
% 此程序为整理面波数据，为克里金插值做准备；
% 输入为读取的面波数据；
% 输出为面波数据点坐标和视横波速度值。
data_sh = data.textdata;  
k = strfind(data_sh,'Location:');
nlie = length(cell2mat(k));
data_sh_length = length(data_sh);

% 数据解译，读出每个频散曲线的起点与长度
%  初始化矩阵
list_begin = ones(1,nlie);
xlocation = ones(1,nlie);
n = 1;
for i = 1:data_sh_length
    if k{i}
        begin = i+1;
        while k{begin}
            begin = begin+1;
        end
        list_begin(n) = begin;
        xlocation(n) = str2double(data_sh{i,2});
        n = n+1;
    end
end

% 创建克里金插值矩阵
points_length = data_sh_length - nlie - 1;
points = zeros(points_length,2);
vs_value = zeros(points_length,1);
nn = 1;
for i = 1:nlie-1
    A = data_sh(list_begin(i):...
        list_begin(i+1)-2,:);
    A = cellcell2mat(A);
    for j = 1:length(A)
        points(nn,1) = xlocation(i);
        points(nn,2) = A(j,1);
        vs_value(nn) = A(j,2);
        nn = nn + 1;
    end
end
clear A
A = data_sh(list_begin(nlie):...
    end,:);
A = cellcell2mat(A);
points(nn:end,1) = xlocation(nlie);
points(nn:end,2) = A(:,1);
vs_value(nn:end) = A(:,2);

3.2 cellcell2mat函数

function C = cellcell2mat(A)
% 此程序为将嵌套元胞数据转为矩阵
[row,col] = size(A);
C = ones(row,col);
for i = 1:col
    for j = 1:row
        a = cell2mat(A(j,i));
        b = str2double(a);
        C(j,i) = b;
    end
end

3.3 sInterp函数

function [xa,ya] = sInterp(xlocation,data_gc,interp_num,num_points)
a = polyfit(xlocation,data_gc,interp_num);
warning('off');
xa = linspace(min(xlocation),max(xlocation),num_points);
ya = polyval(a,xa);

3.4 sLcLabel函数

function data_lclabel = sLcLabel(data_lc)
n = length(data_lc);
data_lclab = num2str(data_lc);
data_lclabel = cell(n,1);
ak1 = data_lclab(1,:);
data_lclabel{1} = strcat(ak1(1:end-3),'+',ak1(end-2:end));
clear ak1 ak2
for i = 2:n
    ak1 = data_lclab(i-1,:);
    ak2 = data_lclab(i,:);
    if strcmp(ak1(1:end-3),ak2(1:end-3))
        data_lclabel{i} = ak2(end-2:end);
    else
        data_lclabel{i} = strcat(ak2(1:end-3),...
        '+',ak2(end-2:end));
    end
end

4、完整代码

close all
clear 
clc

% 此程序功能是面波数据快速出图
% 作者：shangxiang
% 时间：2023年2月23日

% 读取面波数据
[FileName,PathName] = uigetfile('*.txt','请选择视横波速度文件',...
    'MultiSelect','on');
filename = strcat(PathName,FileName);
data = importdata(filename);
fprintf('
读取视横波速度完成!
请按任意键继续...
');

% pause;
% 读取GPS测量高程数据
clear FileName PathName
[FileName,PathName] = uigetfile('*.txt','请选择GPS高程文件',...
    'MultiSelect','on');
filename_gc = strcat(PathName,FileName);
data_gc = load(filename_gc);
fprintf('
读取高程数据完成!
请按任意键继续...
');

% 读取GPS测量里程数据
clear FileName PathName
[FileName,PathName] = uigetfile('*.txt','请选择GPS里程文件',...
    'MultiSelect','on');
filename_lc = strcat(PathName,FileName);
data_lc = load(filename_lc);
fprintf('
读取里程数据完成!
');

% 初始参数设置
% 最大深度
depth_max = 80;
% 插值点数
num_points = 40;

% 面波数据预处理
[points,vs_value,xlocation] = dealMBdata(data);

% 开始克里金插值
% 克里金插值预设参数
theta = [1 1];
lob = [0.1 0.1];
upb = [2 2];
[points_new,vs_value_new] = dsmerge(points,vs_value);
[dmodel,perf] = dacefit(points_new,vs_value_new,@regpoly0,...
    @correxp,theta,lob,upb);
% [dmodel,perf] = dacefit(points,vs_value,@regpoly0,@correxp,theta,lob,upb);
xmin = min(points(:,1));
xmax = max(points(:,1));
XX = gridsamp([xmin 0;xmax depth_max],num_points);
[YX,MSE] = predictor(XX,dmodel);
X1 = reshape(XX(:,1),num_points,num_points);
Y1 = reshape(XX(:,2),num_points,num_points);
YX = reshape(YX,size(X1));

% 对地形数据进行插值，默认插值点数为9，可更改；
interp_num = 9;
[xa,ya] = sInterp(xlocation,data_gc,interp_num,num_points);
Y1_new = -Y1 + ya;

% 对图形进行处理，补充图像下部
% ynew = max(Y1_new(end,:));
% Y1_new = Y1_new(Y1_new > ynew);
% X1 = X1(Y1_new > ynew);
% YX = YX(Y1_new > ynew);
% 处理里程数据
% 获取横坐标位置
data_lcx = data_lc - min(data_lc);
% 获取横坐标刻度
data_lclabel = sLcLabel(data_lc);

% 画图
figure(1);
clear k
k = (depth_max+max(ya))/max(X1(1,:));
set(gcf,'position',[50 150 1200 1500*k]);
% pcolor(X1,Y1_new,YX);
contourf(X1,Y1_new,YX,50,'linecolor','none');
set(gca,'xtick',data_lcx,'xticklabel',...
    data_lclabel,'xticklabelrotation',45);
caxis([100 500]);
colormap(jet);
h = colorbar;
set(get(h,'title'),'string','fontname{宋体}视横波速度（米/秒）',...
    'FontSize',10);
clear a b
axis equal;
box off;
% axis off
shading flat
set(gca,'fontname','times new roman','fontsize',...
    10,'fontweight','normal');
xlabel('fontname{宋体}里程（m）');
ylabel('fontname{宋体}高程（m）');


% 保存数据
clear xx yy zz
xx = X1(:);
yy = Y1_new(:);
zz = YX(:);
C = [xx,yy,zz];
dlmwrite(strcat(PathName,'mianbo.dat'),C);
clear yy
yy = Y1_new(1,:);
high = [xa',yy'];
dlmwrite(strcat(PathName,'gaocheng.dat'),high);

代码运行过程中如果出现bug，请依据实际工程修改。