www.pudn.com > fk_migration.rar > plotseismo.m, change:2015-12-17,size:11006b


 
function [para,datarange,datanew,tnew,xnew]=plotseismo(data,t,x,para); 
%function [para,datarange,datanew,tnew,xnew]=plotseismo(data,t0,dt,x0,dx,para); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%function [para,datarange,datanew,tnew,xnew]=plotseismo(data,t,x,para); 
%功能:用于绘制地震剖面图 
%输入%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%如果无输入,不会绘图,但可以得到para的默认值 
%data...二维数组,地震剖面数据,每一列为一个地震道的数据 
%t0为每一地震道数据对应的时间坐标的起始坐标 
%dt为每一地震道数据对应的时间坐标的间隔 
%x0为各个地震道对应的坐标的起始坐标 
%dx为各个地震道对应的坐标的间隔 
%t...序列,长度同data的行数,每一地震道数据对应的的坐标<默认值:t(0)=0,dt=1,……> 
%x...序列,长度同data的列数,各个地震道对应的坐标<默认值:x(0)=1,dx=1,……> 
%para...详细参数设置 <默认值:para={1,'k',1,'k','w','k',3,0,1.5,0.93,[],1};> 
%       para{1}...linemode(1,显示线;0,不显示)<默认值=1> 
%       para{2}...linecolor线色<默认值='k'(黑色[0 0 0])> 
%       para{3}...shademode 0,无填充; 
%                           1,正单色填充;-1负单色填充;填充色由para(4)决定 
%                           2,正渐变填充;-2负渐变填充;填充色由para(12)决定,此时linecolor取para(12){1} 
%                           其它正,正渐变填充;其它负,负渐变填充;填充色由para(12)决定,此时线色同linecolor 
%       para{4}...shadecolor填充色<默认值='k'(黑色[0 0 0])> 
%       para{5}...bgcolor背景色<默认值='w'(白色[1 1 1])> 
%       para{6}...axescolor<默认值='k'(黑色[0 0 0])> 
%       para{7}...interpo地震道数据插值,每两个原始点之间插值点数<默认值=3> 
%       para{8}...xjump地震道数过多时,越过xjump个道再显示<默认值=0> 
%       para{9}...ampscale,显示时,最大振幅是道间距的ampscale倍<默认值=1.5> 
%       para{10}...shadescale,填充比例<默认值=0.93> 
%       para{11}...datacut,长度为2的序列,原始数据剪切<默认值=[datamin,datamax]> 
%       para{12}...colorrange,长度为2的cell,形式如{[1 1 0.3],[1 0 0]},para(3)~=0,1,-1时起作用, 
%                                                  para(12){1}最小值颜色,para(12){2}最大值颜色 
%                              或者 
%                              一个数,只能是1,2,3,对应三种预设如下 
%                              1:colorrange={[1 1 0.3],[1 0 0]};bgcolor=[0 1 1];axescolor='b'; 
%                              2:colorrange={[0 0 1],[1 1 1]};bgcolor=[0 0 0];axescolor=[0 1 1]; 
%                              3:colorrange={[0 1 1],[0 0 1]};bgcolor=[1 0.9 0.1];axescolor=[1 0 0.2]; 
%                              此预设只在para(3)~=0,1,-1时启用, 
%                              <默认值=1> 
%输出%%%%%%%%%%%%%%%%%%% 
%para...当前绘图所采用的参数 
%datarange...原始数据范围,为下次成图时的datacut提供参考 
%datanew...当前绘图所使用的数据,它是执行了插值和越道后的数据 
%tnew...datanew每一地震道数据对应的的坐标 
%xnew...datanew各个地震道对应的坐标 
% 
%杨昊,波动理论与成像技术实验室,吉林大学,2005 
%本程序仅供吉林大学波动理论与成像技术实验室的成员使用,其他人未经编者允许请勿使用 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%没有输入参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if(nargin==0) 
    disp(['%%%%%%%%%%%%%%%%%%%%%%%%%%']); 
    disp(['没有输入数据,不能绘图。']); 
    disp(['函数返回para的默认值:']); 
    disp(['para{1} = 1  (linemode)']); 
    disp(['para{2} =''k'' (linecolor)']); 
    disp(['para{3} = 1  (shademode)']); 
    disp(['para{4} =''k'' (shadecolor)']); 
    disp(['para{5} =''w'' (bgcolor)']); 
    disp(['para{6} =''k'' (axescolor)']); 
    disp(['para{7} = 3  (interpo)']); 
    disp(['para{8} = 0  (xjump)']); 
    disp(['para{9} =1.5 (ampscale)']); 
    disp(['para{10}=0.93(shadescale)']); 
    disp(['para{11}=[]  (datacut)']); 
    disp(['para{12}= 1  (colorrange)']); 
    disp(['%%%%%%%%%%%%%%%%%%%%%%%%%%']); 
    para={1,'k',1,'k','w','k',3,0,1.5,0.96,[],1}; 
    return; 
%有输入参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
else 
%%%% 
[ndatarow,ndatacol]=size(data); 
datamax=max(max(data)); 
datamin=min(min(data)); 
datarange=[datamin,datamax]; 
if(nargin<2||isempty(t)) 
   t=0:ndatarow-1; 
end 
if(nargin<3||isempty(x)) 
    x=1:ndatacol; 
end 
%if(nargin<3||isempty(t0)) 
 %   t0=0; 
  %  dt=1; 
%end 
%if(nargin<5||isempty(x0)) 
 %   x0=0; 
  %  dx=1; 
%end 
%t=t0+(0:ndatarow-1)*dt; 
%x=x0+(0:ndatacol-1)*dx; 
nt=length(t); 
nx=length(x); 
if(nt~=ndatarow||nx~=ndatacol) 
    error('数据和数据坐标长度不匹配!'); 
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if(nargin<4) 
 %if(nargin<6)   
    para=cell(1,12); 
end 
%%%%%%%%%%%%%%%%% 
linemode=para{1}; 
linecolor=para{2}; 
shademode=para{3}; 
shadecolor=para{4}; 
bgcolor=para{5}; 
axescolor=para{6}; 
interpo=para{7}; 
xjump=para{8}; 
ampscale=para{9}; 
shadescale=para{10}; 
datacut=para{11}; 
colorrange=para{12}; 
%%%%%%%%%%%%%%%%%%%%%% 
if(isempty(linemode)); 
    linemode=1; 
end 
if(isempty(linecolor)); 
    linecolor='k'; 
end 
if(isempty(shademode)); 
    shademode=1; 
end 
if(isempty(shadecolor)); 
    shadecolor='k'; 
end 
if(isempty(bgcolor)); 
    bgcolor='w'; 
end 
if(isempty(axescolor)); 
    axescolor='k'; 
end   
if(isempty(interpo)); 
    interpo=3; 
end 
if(isempty(xjump)); 
    xjump=0; 
end 
if(isempty(ampscale)); 
    ampscale=1.5; 
end 
if(isempty(shadescale)); 
    shadescale=0.93; 
end 
if(isempty(datacut)); 
    datacut=[datamin,datamax]; 
end 
if(isempty(colorrange)); 
    colorrange=1; 
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if(shademode~=0&&abs(shademode)~=1&&length(colorrange)==1) 
    switch(colorrange) 
        case(1) 
            colorrange={[1 1 0.6],[1 0 0]}; 
            bgcolor=[0 1 1]; 
            axescolor='b'; 
        case(2) 
            colorrange={[0 0 1],[1 1 1]}; 
            bgcolor=[0 0 0]; 
            axescolor=[0 1 1]; 
        case(3) 
            colorrange={[0 1 1],[0 0 1]}; 
            bgcolor=[1 0.9 0.1]; 
            axescolor=[1 0 0.2]; 
        otherwise 
            error('para(12)错误,不合法的colorrange值!') 
    end 
end         
if(abs(shademode)==2) 
    linecolor=colorrange{1}; 
end 
%插值和越道%%%%%%%%%%%%%% 
ntnew=(nt-1)*interpo+nt; 
tnew=linspace(t(1),t(end),ntnew); 
nxall=1:(xjump+1):nx; 
xnew=x(nxall); 
datanew=zeros(ntnew,length(xnew)); 
nxnew=0; 
for (ix=nxall) 
   nxnew=nxnew+1; 
   datanew(:,nxnew)=spline(t,data(:,ix),tnew); 
end 
%数据剪切%%%%%% 
nxall=1:nxnew; 
ntall=1:ntnew; 
if(datacut(1)>datamin) 
    itcut=find(datanew<datacut(1)); 
    datanew(itcut)=datacut(1); 
end 
if(datacut(2)<datamax) 
    itcut=find(datanew>datacut(2)); 
    datanew(itcut)=datacut(2); 
end 
%%%%%%%%%%%%%%%%%%%%%% 
dxnew=xnew(2)-xnew(1); 
maxabs=max(max(abs(datanew))); 
plotamp=ampscale*dxnew; 
datanew=plotamp*datanew/maxabs; 
%查找填充域%%%%%%%%% 
sg=sign(shademode); 
if(sg~=0) 
    shadelim=(1-shadescale)*plotamp; 
    nfillrow=0;    
    for(ix=nxall) 
        ablenextrow=1; 
        for(it=ntall)             
            if(ablenextrow)                     
                if(sg*datanew(it,ix)>shadelim) 
                    nfillrow=nfillrow+1; 
                    fillmatrix(nfillrow,1)=ix; 
                    fillmatrix(nfillrow,2)=it; 
                    ablenextrow=0; 
                end 
            else 
                if(sg*datanew(it,ix)<=shadelim) 
                    fillmatrix(nfillrow,3)=it-1; 
                    ablenextrow=1;                         
                end 
            end 
        end 
        if(~ablenextrow) 
            fillmatrix(nfillrow,3)=it; 
        end 
    end 
end 
%关于图的一些基本属性设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
figure1=figure('Color',bgcolor); 
axes1=axes(... 
  'Color',bgcolor,... 
  'XColor',axescolor,... 
  'XAxisLocation','top',... 
  'YColor',axescolor,... 
  'YDir','reverse',... 
  'Parent',figure1); 
axis(axes1,[xnew(1)-plotamp,xnew(end)+plotamp,tnew(1) tnew(end)]); 
hold(axes1,'all'); 
%画线%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if(linemode~=0) 
    for(ix=nxall) 
        plot(datanew(:,ix)+xnew(ix),tnew,'color',linecolor); 
    end 
end 
%填充%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if(sg~=0)     
    for(ifillrow=1:nfillrow) 
        ix=fillmatrix(ifillrow,1); 
        its=fillmatrix(ifillrow,2); 
        itb=fillmatrix(ifillrow,3); 
        %%%%%%%%%%%%%%%%%%%%%%%%%%% 
        lengthfill=3+itb-its; 
        fillx=zeros(1,lengthfill); 
        filly=fillx; 
        sgsh=sg*shadelim; 
        %%%%%%%%%%%%%%%%%%%%%%%%%%% 
        if(its==1) 
            fillx(1)=sgsh; 
            filly(1)=tnew(1); 
        else 
            p1=[tnew(its-1),datanew(its-1,ix)]; 
            p2=[tnew(its),datanew(its,ix)]; 
            fillp=xinterp(p1,p2,sgsh); 
            fillx(1)=fillp(2); 
            filly(1)=fillp(1); 
        end 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        fillx(2:lengthfill-1)=datanew(its:itb,ix); 
        filly(2:lengthfill-1)=tnew(its:itb); 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        if(itb==ntnew) 
            fillx(lengthfill)=sgsh; 
            filly(lengthfill)=tnew(end); 
        else 
            p1=[tnew(itb),datanew(itb,ix)]; 
            p2=[tnew(itb+1),datanew(itb+1,ix)]; 
            fillp=xinterp(p1,p2,sgsh); 
            fillx(lengthfill)=fillp(2); 
            filly(lengthfill)=fillp(1); 
        end 
        %单色填充%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        if(abs(shademode)==1); 
            fill(fillx+xnew(ix),filly,shadecolor,'EdgeColor',linecolor); 
        %渐变色填充%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        else 
            cr1=colorrange{1}; 
            cr2=colorrange{2}; 
            cr2s1=cr2-cr1; 
            pss=plotamp-shadelim; 
            fillx4p=zeros(1,4); 
            filly4p=fillx4p; 
            for(ifi=1:lengthfill-1) 
                p2x=fillx(ifi); 
                p2y=filly(ifi); 
                p3x=fillx(ifi+1); 
                p3y=filly(ifi+1); 
                fillx4p=[sgsh,p2x,p3x,sgsh]; 
                filly4p=[p2y,p2y,p3y,p3y]; 
                crfactor=(abs(p2x+p3x)/2-shadelim)/pss; 
                colorlet=cr1+cr2s1*crfactor; 
                fill(fillx4p+xnew(ix),filly4p,colorlet,'edgecolor',colorlet); 
            end 
        end 
    fillx=[]; 
    filly=[]; 
    end 
end 
%%% 
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
disp(['%%%%%%%%%%%%%%%%%%%%%%%%%%']); 
disp(['plotseismo已经完成绘图工作。']); 
disp(['%%%%%%%%%%%%%%%%%%%%%%%%%%']); 
para={linemode,linecolor,shademode,shadecolor,bgcolor,axescolor,interpo,xjump,ampscale,shadescale,datacut,colorrange}; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%子函数     
function p0=xinterp(p1,p2,y0)%插出与对应y0的x0   
x0=p1(1)+(p2(1)-p1(1))*(y0-p1(2))/(p2(2)-p1(2)); 
p0=[x0,y0]; 
%2005-10-24大功告成,用时两天%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%