chongbianxie.zip_压缩感知 地震_地震 重建_地震压缩感知_地震数据_重建

clc clear all tic num_l=256; %信号长度 for i=1:1:64 beq(:,i)=rick1(15,0.25)+0.8*rick1(18,0.5+0.0001*i*i); %合成地震记录 end figure(1) wigb(beq) xlabel("trace") ylabel("time") qs_num=20;%缺失道 %qsd=randi(64,1,20);%缺失道号 qsd=4:3:61; %qsd=sort(qsd,'ascend'); beq1=beq; beq2=beq; beq1(:,qsd)=zeros(num_l,qs_num); figure(2) wigb(beq1) beq3=beq1; xlabel("trace") ylabel("time") caiy=ones(num_l,64); caiy(:,qsd)=zeros(num_l,qs_num); figure(3) wigb(caiy) xlabel("trace") ylabel("time") beq1(:,qsd)=[]; figure(4) wigb(beq1) xlabel("trace") ylabel("time") hh=1; for i=1:1:32 for j=1:1:5 yy(:,hh)=reshape(beq1((i-1)*8+1:i*8,(j-1)*8+1:j*8),1,[]); hh=hh+1; end end %%%%%%%%%%%%%%%%%%%%%c产生初始字典 K=256; Pn=ceil(sqrt(K)); bb=8; DCT=zeros(bb,Pn); for k=0:1:Pn-1 V=cos([0:1:bb-1]'*k*pi/Pn); if k>0,V=V-mean(V); end DCT(:,k+1)=V/norm(V); end G=kron(DCT,DCT); %%%%%%%%%%%%%%%%显示字典 for i=1:1:K M(1:8,(i-1)*8+1:i*8)=reshape(G(:,i),8,[]); end k=1; for i=1:1:sqrt(K) for j=1:1:sqrt(K) N((i-1)*8+1:i*8,(j-1)*8+1:j*8)=M(:,(k-1)*8+1:k*8); k=k+1; end end figure(5) wigb(N) xlabel("trace") ylabel("time") size(N); %%%%%%%%%%%%%%%%%%%%%%%%%%%开始字典学习 yb_nu=128; %所取样本个数 c_max=4; %---------选取重构原子的个数--- l_nu=100; %--------- -设置学习的次数 for num=1:1:l_nu for h=1:1:yb_nu x1=zeros(K,yb_nu); y=yy(:,h); for i=1:c_max for m=1:K inner(m)=dot(y,G(:,m)); end [B3,IX3]=sort(inner,'descend');%降序排列,B3为排序后 %的向量，IX3为索引 x1(IX3(1),h)=(B3(1))/(norm(G(:,IX3(1))))^2; %向量的2范数，即长度 y=y-x1(IX3(1),h)*G(:,IX3(1)); %计算残差 end end for i=1:1:K for j=1:1:yb_nu jj=1; if x1(i,j)~=0 im(jj)=j; jj=jj+1; end end if jj==1 break end w=zeros(yb_nu,jj-1); for k=1:1:jj-1 w(im(k),k)=1; end yy1=yy*w; x1_k=x1*w;%去除系数中的零向量 Ek=yy1-G*x1_k+G(:,i)*x1_k(i,:); [U B V]=svd(Ek); %对E进行scd分解 G(:,i)=U(:,1); %更新字典 x1(i,im)=V(:,1)*B(1,1); %更新该原子对应的系数 end end %%%%%%%%%%%%%%%%结束字典学习过程 fid=fopen('学习后的字典.txt','w'); HSA1=size(G,1); fprintf(fid,'%s'); fprintf(fid,'\r\n'); for i=1:HSA1 fprintf(fid,'%10.2f',G(i,:)); fprintf(fid,'\r\n'); end G=load('学习后的字典.txt'); %%%%%%%%%%%%%%重构过程 mm=1; X=zeros(256,256); for i=1:1:32 for j=1:1:8 kk=1; km=[]; for k=1:1:8 if caiy((i-1)*8+1,(j-1)*8+k)==1 km(kk)=k; kk=kk+1; end end cs=zeros(8*(kk-1),64); a=eye(8); for m=1:1:kk-1 cs((m-1)*8+1:8*m,(km(m)-1)*8+1:8*km(m))=a; % 产生采样矩阵 end gc=zeros((kk-1)*8,1); %观测值初始化 gc=cs*reshape(beq((i-1)*8+1:i*8,(j-1)*8+1:j*8),1,[])'; %生成观测值 ysgz=cs*G; %压缩传感矩阵 X(:,mm)=StOMP(gc,ysgz,c_max ); %运用StROMP计算分解系数 mm=mm+1; end end size(X) beq_cg=G*X; for i=1:1:256 MM(1:8,(i-1)*8+1:i*8)=reshape(beq_cg(:,i),8,[]); end k=1; for i=1:1:32 for j=1:1:8 NN((i-1)*8+1:i*8,(j-1)*8+1:j*8)=MM(:,(k-1)*8+1:k*8); k=k+1; end end figure(6) wigb(NN); size(NN) figure(7) error=beq-NN; wigb(error); min(min(error)); max(max(error)); %------ ----计算误差------ --- error_1=reshape(error,1,[]); beq_1=reshape(beq,1,[]); beq_qs=reshape(beq3,1,[]); NN1=reshape(NN,1,[]); mse1=error_1*error_1'/(256*64) %原始与重构的均方差 error_2=beq_1-beq_qs; snr1=10*log10(beq_1*beq_1'/(error_2*error_2')) %原始与缺失的信噪比 error_3=beq_1-NN1; snr1=10*log10(beq_1*beq_1'/(error_3*error_3')) %原始与重构的信噪比