程序我运行了没有问题，但是画的图上面有个vpR的标志，我在程序里没有看到这个语句，不知道为什么会出现，希望能去掉。

clear;

clc;

close all;

vp1=1000/171;

vs1=1000/256;

vp2=1000/230;

vs2=1000/394;

vp3=1000/617;

vs3=0;

ro1=2.65;

ro2=2.6;

ro3=1.05;

vqu=0.8;

vsh=0.2;

kqu=ro1*vp1*vp1-4*ro1*vs1*vs1/3;

gqu=ro1*vs1*vs1;

ksh=ro2*vp2*vp2-4*ro2*vs2*vs2/3;

gsh=ro2*vs2*vs2;

kup=kqu*vqu+ksh*vsh;

klw=1./(vqu/kqu+vsh/ksh);

gup=gqu*vqu+gsh*vsh;

glw=1./(vqu/gqu+vsh/gsh);

ko=(kup+klw)/2;

go=(gup+glw)/2;

kfl=ro3*vp3*vp3;

figure

i=1;

for pt=0.01:0.01:0.4

a1=0.01;

ro_p=((1-pt)*ro1+(1-pt)*ro2)/2+pt*ro3;

[Tsh,Fsh]=XuWhitePQ(ko,go,0,0,a1);

[Tqu,Fqu]=XuWhitePQ(ko,go,0,0,a1);

p=(vsh.*Tsh+vqu.*Tqu)/3;

q=(vsh.*Fsh+vqu.*Fqu)/5;

kdry=ko.*(1-pt).^p;

gdry=go.*(1-pt).^q;

ksat=kdry+(1-kdry./ko).^2 ./ ( pt./kfl + (1-pt)./ko + kdry./ko./ko +1.0e-12);

gsat=gdry;

vpsat(i,1)=sqrt((ksat+4*gsat/3)*1.0e6./ro_p);

vssat(i,1)=sqrt( gsat*1.0e6./ro_p );

pt1(i,1)=pt;

i=i+1;

end

plot(pt1,vpsat,'-rs','LineWidth',1.5);hold on;

plot(pt1,vssat,'-gs','LineWidth',1.5);hold on;

% h1 = legend('vp','vs',2);  % 2 为 NorthWest 请help legend

% set(h1,'Interpreter','none')

text(pt1(1,1),vpsat(1,1),'0.01','FontSize',16)

text(pt1(1,1),vssat(1,1),'0.01','FontSize',16)

i=1;

for pt=0.01:0.01:0.4

a1=0.05;

ro_p=((1-pt)*ro1+(1-pt)*ro2)/2+pt*ro3;

[Tsh,Fsh]=XuWhitePQ(ko,go,0,0,a1);

[Tqu,Fqu]=XuWhitePQ(ko,go,0,0,a1);

p=(vsh.*Tsh+vqu.*Tqu)/3;

q=(vsh.*Fsh+vqu.*Fqu)/5;

kdry=ko.*(1-pt).^p;

gdry=go.*(1-pt).^q;

ksat=kdry+(1-kdry./ko).^2 ./ ( pt./kfl + (1-pt)./ko + kdry./ko./ko +1.0e-12);

gsat=gdry;

vpsat(i,1)=sqrt((ksat+4*gsat/3)*1.0e6./ro_p);

vssat(i,1)=sqrt( gsat*1.0e6./ro_p );

pt1(i,1)=pt;

i=i+1;

end

plot(pt1,vpsat,'-rp','LineWidth',1.5);hold on;

text(0.2,4000,'vp\rightarrow')

plot(pt1,vssat,'-gp','LineWidth',1.5);hold on;

text(pt1(10,1),vpsat(10,1),'0.05','FontSize',16)

text(pt1(40,1),vssat(40,1),'0.05','FontSize',16)

% h2 = legend('vp','vs',2);  % 2 为 NorthWest 请help legend

% set(h2,'Interpreter','none')

i=1;

for pt=0.01:0.01:0.4

a1=0.1;

ro_p=((1-pt)*ro1+(1-pt)*ro2)/2+pt*ro3;

[Tsh,Fsh]=XuWhitePQ(ko,go,0,0,a1);

[Tqu,Fqu]=XuWhitePQ(ko,go,0,0,a1);

p=(vsh.*Tsh+vqu.*Tqu)/3;

q=(vsh.*Fsh+vqu.*Fqu)/5;

kdry=ko.*(1-pt).^p;

gdry=go.*(1-pt).^q;

ksat=kdry+(1-kdry./ko).^2 ./ ( pt./kfl + (1-pt)./ko + kdry./ko./ko +1.0e-12);

gsat=gdry;

vpsat(i,1)=sqrt((ksat+4*gsat/3)*1.0e6./ro_p);

vssat(i,1)=sqrt( gsat*1.0e6./ro_p );

pt1(i,1)=pt;

i=i+1;

end

plot(pt1,vpsat,'-rd','LineWidth',1.5);hold on;

plot(pt1,vssat,'-gd','LineWidth',1.5);hold on;

text(pt1(36,1),vpsat(36,1),'0.1','FontSize',16)

text(pt1(40,1),vssat(40,1),'0.1','FontSize',16)

% h3 = legend('vp','vs',2);  % 2 为 NorthWest 请help legend

% set(h3,'Interpreter','none')

i=1;

for pt=0.01:0.01:0.4

a1=0.15;

ro_p=((1-pt)*ro1+(1-pt)*ro2)/2+pt*ro3;

[Tsh,Fsh]=XuWhitePQ(ko,go,0,0,a1);

[Tqu,Fqu]=XuWhitePQ(ko,go,0,0,a1);

p=(vsh.*Tsh+vqu.*Tqu)/3;

q=(vsh.*Fsh+vqu.*Fqu)/5;

kdry=ko.*(1-pt).^p;

gdry=go.*(1-pt).^q;

ksat=kdry+(1-kdry./ko).^2 ./ ( pt./kfl + (1-pt)./ko + kdry./ko./ko +1.0e-12);

gsat=gdry;

vpsat(i,1)=sqrt((ksat+4*gsat/3)*1.0e6./ro_p);

vssat(i,1)=sqrt( gsat*1.0e6./ro_p );

pt1(i,1)=pt;

i=i+1;

end

plot(pt1,vpsat,'-rh','LineWidth',1.5);hold on;

text(pt1(20,1),vpsat(20,1),'0.15','FontSize',16)

plot(pt1,vssat,'-gh','LineWidth',1.5);hold on;

text(pt1(20,1),vssat(20,1),'0.15','FontSize',16)

text(0.05,1000,'vs \rightarrow')

% h4= legend('vp','vs',2);  % 2 为 NorthWest 请help legend

% set(h4,'Interpreter','none')

% title('纵横波速度与孔隙度关系')

xlabel('孔隙度');

ylabel('速度 (m/s)');

111.png

(1.62 KB, 下载次数: 0)

2014-5-21 15:04 上传

图上的vpr我不知道是怎么出来的