三维k-means无监督智能地震速度拾取（matlab实现）

本文从三个方向探究了叠加速度拾取智能化的方法。
首先是速度谱预处理方面，经过统一阈值筛选后的速度谱往往忽略了浅层的能量团，而且深层速度的聚焦性差。故设计了一种随深度的变化，阈值也相应变化的滑动窗来实现阈值的灵活筛选，使浅层的速度也能被拾取，深层速度的拾取更加准确。
其次是聚类算法的升级，传统的二维k means聚类算法，得的中心是截面的几何中心,而非真正速度谱的能量中心。特创新出三维k means聚类，使拾取的结果更加接近于真实的能量团，从而得自动更加准确的叠加速度来进行动校正。
最后是其他聚类算法对聚类结果的比较，k means聚类算法需要人工设置聚类中心个数和最大循环次数，db-scan聚类算法不需要设置聚类中心个数，只需要调节限制圆的半径和圆内最小点数，在算法运行的过程中就可以获得聚类中心个数。在本次科创过程中探究了这两种算法的优缺点更好地改进对速度拾取算提出更好的改进措施。

叠加速度拾取原理

利用速度谱方法求取叠加速度是目前生产单位提取速度参数的重要手段。所谓叠加速度，是指对道集内某个反射波同相轴用不同的速度进行动校正并分析校正后的叠加效果，其中叠加效果最好的那个速度就是该反射波的叠加速度。
为了得到叠加速度，需要将原始地震道集和界面反射系数褶积得到叠加速度谱，用阈值对速度谱进行筛选，用算法对筛选后的速度谱进行聚类处理，处理得到的聚类中心的坐标对应的就是不同深度的叠加速度大小。
用拾取的地震速度对地震道集进行动校正，就可以将双曲线拉直，再进行后期的水平叠加。

K-means聚类原理简介

1.输入一聚类簇数为k的样本集，样本集大小为m*n的矩阵在样本集中随机抽取k个样本作为初始均值向量（cluster center）；
2.计算每个样本点与各均值向量的距离d 根据距离d将每个样本点划入与之最近的簇中，此时均值向量发生改变；
3.对簇中每个样本求和/簇中样本个数求得新的均值向量；
4.更新均值向量，当均值向量不再改变时停止迭代，输出将样本集最后划分的k个簇；
为了探究出一组最合适的参数，我们调节了不同的阈值和聚类中心个数，研究其对聚类效果的影响:
从上图可以明显的发现，阈值过大会丢失大部分的能量团数据使浅层的速度拾取效果不佳；阈值过不断地调试大量的噪音。经过不断的调试最后选择95为最佳阈值。

从上图可以发现不同聚类中心的个数也会影响聚类效果，其中12个聚类中心的效果是最佳的。

自适应滑窗预处理

用统一阈值对速度谱进行筛选时固定的阈值使其拾取精度在一定程度上受到速度能量团聚焦性的制约，而且容易忽略弱反射速度的能量团。故在筛选速度谱时，加入随时间变化的滑动窗，来实现对能量的灵活筛选。我们设计的滑动窗口的窗函数为
（其中i为第i时刻，d-win为窗长）。随着时间i的增大，滑动窗口从速度谱的浅部向深部移动，每一时刻窗内的阈值为
。（其中都为窗长）。可以发现滑动窗内的阈值会随时刻i的变化而变化，这就实现了阈值的灵活改变。

从图中可以看出加了自适应滑窗后的速度谱的识别效果在浅层和深层明显更优于无滑窗预处理的拾取效果。

二/三维聚类效果对比

传统的k-means算法拾取速度谱是对速度谱设定统一的阈值，大于阈值的部分设定为1，小于阈值的部分设定为0。再对筛选后的二维速度谱进行聚类，由于聚类的对象是二维平面，所以聚类的标准——距离的公式也是二维的距离公式即
二维k-means聚类后的结果是速度谱以阈值做截面后的几何中心而不是真正的能量中心。创新后的三维k-means算法将经过统一阈值筛选后的0,1速度谱再乘以原来的速度谱得到三维的速度谱，再该谱进行聚类。聚类过程中的距离公式变为

三维聚类后的结果在物理意义上相较于二维聚类更加接近于速度谱的能量中心，也更加准确。
结果图如下：
三维k-means聚类算法：

% clear all
% clc
% close all
FONTSIZE=24;   FONTNAME='Times New Roman'; LINEWIDTH=2;
%%%%%%%%%%%%Auto Semblance Spectrum Picking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%forward model
%% 生成速度模型
load stack_velspec.mat; %改变速度谱    速度谱私信我拿
nt=749;
CMP_num=1;
[nt1,n_vel]=size(VEL_SPEC);
x=200:50:2000;
%(1) Thresholding the velocity spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:CMP_num
vel_spec=VEL_SPEC;
vel_spec2=zeros(size(vel_spec));
per_max=95;
per_min=95;
per_interval=per_max-per_min;
d_win=40;
per_dv=200;
start=40;
ultimate=nt1;
thres=97;%%阈值百分数。%求整体阈值的百分比
thres_value = prctile(prctile(vel_spec,thres),thres);  % Set the threshold value 表示被调群体中有n%的数据小于此数值（thres_value是阈值的的具体数）
VEL_SPEC2=(vel_spec>thres_value);                % Filter out small value%(2) Build the training set.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
central_initial=cell(CMP_num,1);
dv=(vmax-vmin)/n_vel;
tic%%%计时
for m=1:CMP_numvel_spec2=VEL_SPEC2.*VEL_SPEC;
vel_spec=VEL_SPEC;
n_point=sum(sum(vel_spec2~=0));
train_set=zeros(n_point,2); %%train-set是不为零的数的位置。 Each example in the training set has two features,i_point=1;                  % Build the training set
for ix=1:n_velfor iz=1:nt1if(vel_spec2(iz,ix)>0)train_set(i_point,1)=ix;train_set(i_point,2)=iz;i_point=i_point+1;endend
end
train_set(all(train_set==0,2),:)=[]
n_point=size(train_set,1)
%(3) Training...
%%%%%%%%%%%%%%%%%%(3.1) Initialize
%%%%%%%%%%%%%%%%%
num_cluster=12;                    % Set the initial number of cluster
dx=floor(n_vel/(num_cluster+1));
dz=floor(nt1/(num_cluster+1));
central=zeros(num_cluster,2);
maxiter=10; %%%%%%%最大循环次数
Max_Stretch=80;for ik=1:num_cluster               % Initialize each cluster centralcentral(ik,1)=dx*ik;central(ik,2)=dz*ik;
end
central_initial{m}=central;%%%%%%%找出开始的聚类中心
% figure;
% imagesc(vel_spec);
% hold on
% plot(central(:,1),central(:,2),'r*');
%(3.2) Interation start...
%%%%%%%%%%%%%%%%%%%%%%%%%%for iter=1:maxiter% Plot each cluster central on the velocity spectrum%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     figure;
%     imagesc(vel_spec)
%     hold on
%     plot(central(:,1),central(:,2),'r*');pause(1);hold off%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%res=0.0;Y=zeros(n_point,1);for ip=1:n_point        % Compute the distance of each points to the first cluster centraldist1=(train_set(ip,1)-central(1,1))^2+(train_set(ip,2)-central(1,2))^2+(vel_spec2(train_set(ip,1),train_set(ip,2))-vel_spec2(central(1,1),central(1,2)))^2;Y(ip)=1;% Compute the distance of each points to the each cluster centralfor ik=2:num_cluster  % Loop over each cluster central %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(特别注意！！！！！！自己改的)   dist_temp=(train_set(ip,1)-central(ik,1))^2+(train_set(ip,2)-central(ik,2))^2+(vel_spec2(train_set(ip,1),train_set(ip,2))-vel_spec2(central(ik,1),central(ik,2)))^2;if(dist1>dist_temp)dist1=dist_temp;Y(ip)=ik;%%%%%%%%%Y矩阵里是数据集不为零的点到每个新中心点的最近距离划分的类别endendres=res+dist1;        % Compute residualend fprintf('iter= %d res = %f\n',iter,res)RES(iter)=res;% Compute the new central of each cluster%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%num_K=zeros(num_cluster,1);central_temp=zeros(size(central));for ip=1:n_pointcentral_temp(Y(ip),1)=central_temp(Y(ip),1)+train_set(ip,1);central_temp(Y(ip),2)=central_temp(Y(ip),2)+train_set(ip,2);num_K(Y(ip))=num_K(Y(ip))+1;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算新的聚类中心for ik=1:num_clusterif(num_K(ik)==0) num_K(ik)=1; endend% Compute the average%%%%%%%%%%%%%%%%%%%%%for ik=1:num_clustercentral_temp(ik,1)=floor(central_temp(ik,1)/num_K(ik));central_temp(ik,2)=floor(central_temp(ik,2)/num_K(ik));end central=central_temp;%%%%%%central 是新的聚类中心central(all(central==0,2),:)=[]num_cluster=size(central,1)central2=centralend%(4) postprocessing% velocity field vnmo=vmin+central2(:,1)*dv;
tnmo=central2(:,2)*dt*dR;
CENTRAL=central2;end
end
toc%%%%%Plot Figure
xx=vmin:dv:vmax;
zz=(1:nt)*dt;figure;
subplot(1,2,1);imagesc(xx,zz,VEL_SPEC2);
xlabel('Velocity(m/s)')
ylabel('Time (s)');ylim([0 1.5]);
set(gca,'tickdir','out','Fontname',FONTNAME,'Fontsize',FONTSIZE);subplot(1,2,2);imagesc(xx,zz,VEL_SPEC);
hold on;
plot(vmin+CENTRAL(:,1)*dv,CENTRAL(:,2)*dt*dR,'r','linewidth',2)
hold on;plot(vmin+CENTRAL(:,1)*dv,CENTRAL(:,2)*dt*dR,'r*');
%plot(vrms(:,K),t0,'g','linewidth',2);
xlabel('Velocity(m/s)')
ylabel('Time (s)')
ylim([0 1.5]);
set(gca,'tickdir','out','Fontname',FONTNAME,'Fontsize',FONTSIZE);`

无监督智能地震速度拾取（matlab实现）相关推荐

无监督构建词库：更快更好的新词发现算法
作者丨苏剑林单位丨追一科技研究方向丨NLP,神经网络个人主页丨kexue.fm 新词发现是 NLP 的基础任务之一,主要是希望通过无监督发掘一些语言特征(主要是统计特征),来判断一批语料中哪些字 ...
从无监督构建词库看「最小熵原理」，套路是如何炼成的
作者丨苏剑林单位丨广州火焰信息科技有限公司研究方向丨NLP,神经网络个人主页丨kexue.fm 在深度学习等端到端方案已经逐步席卷 NLP 的今天,你是否还愿意去思考自然语言背后的基本原理?我们 ...
Google “推翻”无监督研究成果！斩获 ICML 2019 最佳论文
作者 | 夕颜.Just 出品 | AI科技大本营(ID:rgznai100) 6 月 11 日,在美国加州长滩举行的 ICML 公布了 2019 年最佳论文奖,来自苏黎世联邦理工大学.谷歌大脑等的团 ...
基于模型与不基于模型的深度增强学习_CVPR2018: 基于时空模型无监督迁移学习的行人重识别...
Unsupervised Cross-dataset Person Re-identification by Transfer Learning of Spatial-Temporal Pattern ...
服务于离群点检测的无监督特征选择值-特征层次耦合模型
Unsupervised Feature Selection for Outlier Detection by Modelling Hierarchical Value-Feature Couplin ...
【CVPR智慧城市挑战赛】无监督交通异常检测，冠军团队技术分享
[新智元导读]"智能交通视频分析界的ImageNet竞赛"--英伟达城市挑战赛落下帷幕.新加坡松下研究院联合中科院自动化所,提出了一种双模态动静联合检测方案,在交通异常检测比赛中拔 ...
David P.Williams论文系列基于间隙度的声呐图像快速无监督海底特征描述
摘要提出了一种基于侧扫声呐图像的海底无监督特征提取方法.该方法基于间隙度,通过传感器数据测量像素强度变化.不需要训练数据,不需要对像素的统计分布作任何假设,也不需要列举或知道(离散的)海底类型的背景 ...
9种有监督与3种无监督机器学习算法
机器学习作为目前的热点技术广泛运用于数据分析领域,其理论和方法用于解决工程应用的复杂问题.然而在机器学习领域,没有算法能完美地解决所有问题(数据集的规模与结构.性能与便利度.可解释性等不可能三角),识 ...
区块链小组作业：无界智能运动竞技类APP
无界智能运动竞技类APP 1. 前言 1.1 编写目的 1.2 背景 2. 摘要 3. 总体设计 3.1 需求规定 3.2 "区块链+证书" 3.2.1 业务流程 3.2.2 证书 ...

无监督智能地震速度拾取（matlab实现）