三維k-means無(wú)監(jiān)督智能地震速度拾取（matlab實(shí)現(xiàn)）

本文從三個(gè)方向探究了疊加速度拾取智能化的方法。
首先是速度譜預(yù)處理方面，經(jīng)過(guò)統(tǒng)一閾值篩選后的速度譜往往忽略了淺層的能量團(tuán)，而且深層速度的聚焦性差。故設(shè)計(jì)了一種隨深度的變化，閾值也相應(yīng)變化的滑動(dòng)窗來(lái)實(shí)現(xiàn)閾值的靈活篩選，使淺層的速度也能被拾取，深層速度的拾取更加準(zhǔn)確。
其次是聚類(lèi)算法的升級(jí)，傳統(tǒng)的二維k means聚類(lèi)算法，得的中心是截面的幾何中心,而非真正速度譜的能量中心。特創(chuàng)新出三維k means聚類(lèi)，使拾取的結(jié)果更加接近于真實(shí)的能量團(tuán)，從而得自動(dòng)更加準(zhǔn)確的疊加速度來(lái)進(jìn)行動(dòng)校正。
最后是其他聚類(lèi)算法對(duì)聚類(lèi)結(jié)果的比較，k means聚類(lèi)算法需要人工設(shè)置聚類(lèi)中心個(gè)數(shù)和最大循環(huán)次數(shù)，db-scan聚類(lèi)算法不需要設(shè)置聚類(lèi)中心個(gè)數(shù)，只需要調(diào)節(jié)限制圓的半徑和圓內(nèi)最小點(diǎn)數(shù)，在算法運(yùn)行的過(guò)程中就可以獲得聚類(lèi)中心個(gè)數(shù)。在本次科創(chuàng)過(guò)程中探究了這兩種算法的優(yōu)缺點(diǎn)更好地改進(jìn)對(duì)速度拾取算提出更好的改進(jìn)措施。

疊加速度拾取原理

利用速度譜方法求取疊加速度是目前生產(chǎn)單位提取速度參數(shù)的重要手段。所謂疊加速度，是指對(duì)道集內(nèi)某個(gè)反射波同相軸用不同的速度進(jìn)行動(dòng)校正并分析校正后的疊加效果，其中疊加效果最好的那個(gè)速度就是該反射波的疊加速度。
為了得到疊加速度，需要將原始地震道集和界面反射系數(shù)褶積得到疊加速度譜，用閾值對(duì)速度譜進(jìn)行篩選，用算法對(duì)篩選后的速度譜進(jìn)行聚類(lèi)處理，處理得到的聚類(lèi)中心的坐標(biāo)對(duì)應(yīng)的就是不同深度的疊加速度大小。
用拾取的地震速度對(duì)地震道集進(jìn)行動(dòng)校正，就可以將雙曲線拉直，再進(jìn)行后期的水平疊加。

K-means聚類(lèi)原理簡(jiǎn)介

1.輸入一聚類(lèi)簇?cái)?shù)為k的樣本集，樣本集大小為m*n的矩陣 在樣本集中隨機(jī)抽取k個(gè)樣本作為初始均值向量（cluster center）；
2.計(jì)算每個(gè)樣本點(diǎn)與各均值向量的距離d 根據(jù)距離d將每個(gè)樣本點(diǎn)劃入與之最近的簇中，此時(shí)均值向量發(fā)生改變 ；
3.對(duì)簇中每個(gè)樣本求和/簇中樣本個(gè)數(shù)求得新的均值向量；
4.更新均值向量，當(dāng)均值向量不再改變時(shí)停止迭代 ，輸出將樣本集最后劃分的k個(gè)簇 ；
為了探究出一組最合適的參數(shù)，我們調(diào)節(jié)了不同的閾值和聚類(lèi)中心個(gè)數(shù)，研究其對(duì)聚類(lèi)效果的影響:
從上圖可以明顯的發(fā)現(xiàn)，閾值過(guò)大會(huì)丟失大部分的能量團(tuán)數(shù)據(jù)使淺層的速度拾取效果不佳；閾值過(guò)不斷地調(diào)試大量的噪音。經(jīng)過(guò)不斷的調(diào)試最后選擇95為最佳閾值。

從上圖可以發(fā)現(xiàn)不同聚類(lèi)中心的個(gè)數(shù)也會(huì)影響聚類(lèi)效果，其中12個(gè)聚類(lèi)中心的效果是最佳的。

自適應(yīng)滑窗預(yù)處理

用統(tǒng)一閾值對(duì)速度譜進(jìn)行篩選時(shí)固定的閾值使其拾取精度在一定程度上受到速度能量團(tuán)聚焦性的制約，而且容易忽略弱反射速度的能量團(tuán)。故在篩選速度譜時(shí)，加入隨時(shí)間變化的滑動(dòng)窗，來(lái)實(shí)現(xiàn)對(duì)能量的靈活篩選。我們?cè)O(shè)計(jì)的滑動(dòng)窗口的窗函數(shù)為
（其中i為第i時(shí)刻，d-win為窗長(zhǎng)）。隨著時(shí)間i的增大，滑動(dòng)窗口從速度譜的淺部向深部移動(dòng)，每一時(shí)刻窗內(nèi)的閾值為
。（其中都為窗長(zhǎng)）。可以發(fā)現(xiàn)滑動(dòng)窗內(nèi)的閾值會(huì)隨時(shí)刻i的變化而變化，這就實(shí)現(xiàn)了閾值的靈活改變。

從圖中可以看出加了自適應(yīng)滑窗后的速度譜的識(shí)別效果在淺層和深層明顯更優(yōu)于無(wú)滑窗預(yù)處理的拾取效果。

二/三維聚類(lèi)效果對(duì)比

傳統(tǒng)的k-means算法拾取速度譜是對(duì)速度譜設(shè)定統(tǒng)一的閾值，大于閾值的部分設(shè)定為1，小于閾值的部分設(shè)定為0。再對(duì)篩選后的二維速度譜進(jìn)行聚類(lèi)，由于聚類(lèi)的對(duì)象是二維平面，所以聚類(lèi)的標(biāo)準(zhǔn)——距離的公式也是二維的距離公式即
二維k-means聚類(lèi)后的結(jié)果是速度譜以閾值做截面后的幾何中心而不是真正的能量中心。創(chuàng)新后的三維k-means算法將經(jīng)過(guò)統(tǒng)一閾值篩選后的0,1速度譜再乘以原來(lái)的速度譜得到三維的速度譜，再該譜進(jìn)行聚類(lèi)。聚類(lèi)過(guò)程中的距離公式變?yōu)?/p>

三維聚類(lèi)后的結(jié)果在物理意義上相較于二維聚類(lèi)更加接近于速度譜的能量中心，也更加準(zhǔn)確。
結(jié)果圖如下：
三維k-means聚類(lèi)算法 ：

% 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;%%閾值百分?jǐn)?shù)。%求整體閾值的百分比 thres_value = prctile(prctile(vel_spec,thres),thres); % Set the threshold value 表示被調(diào)群體中有n%的數(shù)據(jù)小于此數(shù)值（thres_value是閾值的的具體數(shù)） 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%%%計(jì)時(shí) 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是不為零的數(shù)的位置。 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; %%%%%%%最大循環(huán)次數(shù) 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;%%%%%%%找出開(kāi)始的聚類(lèi)中心 % 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矩陣?yán)锸菙?shù)據(jù)集不為零的點(diǎn)到每個(gè)新中心點(diǎn)的最近距離劃分的類(lèi)別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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%計(jì)算新的聚類(lèi)中心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 是新的聚類(lèi)中心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);`

總結(jié)

以上是生活随笔為你收集整理的无监督智能地震速度拾取（matlab实现）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。