国家关于网站信息建设管理文件,安陆网站,如何推广视频号,阳江网络问政平台 阳东中学代码分享#xff1a;面波数据快速成图
前言
目前#xff0c;物探数据主要用surfer软件成图#xff0c;surfer软件具有强大的插值和绘图功能#xff0c;成图比较美观。但是#xff0c;生产过程中大量的物探数据#xff0c;依靠excel和surfer来成图耗费人力时间成本。本博…代码分享面波数据快速成图
前言
目前物探数据主要用surfer软件成图surfer软件具有强大的插值和绘图功能成图比较美观。但是生产过程中大量的物探数据依靠excel和surfer来成图耗费人力时间成本。本博文在MATLAB平台上开发了一套用于面波数据快速成图的小程序仅供大家借鉴。
文章目录代码分享面波数据快速成图前言1、成图效果展示1.1 原始图像1.2 高程转换1.3 里程换算1.4 图像加工2、数据读取与图像保存2.1 读取面波视横波速度数据2.2 数据与图像保存3、自编函数3.1 dealMBdata函数3.2 cellcell2mat函数3.3 sInterp函数3.4 sLcLabel函数4、完整代码代码运行过程中如果出现bug请依据实际工程修改。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(\n读取视横波速度完成!\n请按任意键继续...\n);提取数据自编函数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_lengthif k{i}begin i1;while k{begin}begin begin1;endlist_begin(n) begin;xlocation(n) str2double(data_sh{i,2});n n1;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-1A data_sh(list_begin(i):...list_begin(i1)-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:colfor j 1:rowa cell2mat(A(j,i));b str2double(a);C(j,i) b;end
end3.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:nak1 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);elsedata_lclabel{i} strcat(ak2(1:end-3),...,ak2(end-2:end));end
end4、完整代码
close all
clear
clc% 此程序功能是面波数据快速出图
% 作者shangxiang
% 时间2023年2月23日% 读取面波数据
[FileName,PathName] uigetfile(*.txt,请选择视横波速度文件,...MultiSelect,on);
filename strcat(PathName,FileName);
data importdata(filename);
fprintf(\n读取视横波速度完成!\n请按任意键继续...\n);% pause;
% 读取GPS测量高程数据
clear FileName PathName
[FileName,PathName] uigetfile(*.txt,请选择GPS高程文件,...MultiSelect,on);
filename_gc strcat(PathName,FileName);
data_gc load(filename_gc);
fprintf(\n读取高程数据完成!\n请按任意键继续...\n);% 读取GPS测量里程数据
clear FileName PathName
[FileName,PathName] uigetfile(*.txt,请选择GPS里程文件,...MultiSelect,on);
filename_lc strcat(PathName,FileName);
data_lc load(filename_lc);
fprintf(\n读取里程数据完成!\n);% 初始参数设置
% 最大深度
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_maxmax(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请依据实际工程修改。
