二维k-means聚类后的结果是速度谱以阈值做截面后的几何中心而不是真正的能量中心 。创新后的三维k-means算法将经过统一阈值筛选后的0,1速度谱再乘以原来的速度谱得到三维的速度谱，再该谱进行聚类 。聚类过程中的距离公式变为
三维聚类后的结果在物理意义上相较于二维聚类更加接近于速度谱的能量中心，也更加准确 。
结果图如下：
三维k-means聚类算法 ：
【matlab实现无监督智能地震速度拾取】% clear all% clc% close allFONTSIZE=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_numvel_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 = http://www.kingceram.com/post/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;endendendtrain_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 clusterdx=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;endcentral_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;endendtoc%%%%%Plot Figurexx=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);`