biot 双向介质正演模拟matlab程序

biot 双向介质正演模拟matlab程序

function [vp1freqvp2vsq1invq2invqsinv]=biot（vpdryvsdryk0mu0ro0roflkflnuporpermaalfad1d2opt）
%[VP1FREQVP2VSQ1INVQ2INVQSINV]
%      =BIOT（VPDRYVSDRYK0MU0RO0ROFLKFLNUPORPERMAALFAD1D2‘PLOTOPT‘）
%
%	calculates the complete Biot velocity dispersion and attenuation
%	curves for all frequencies. 
% Inputs: VPDRY: P-wave velocity of dry porous rock
%	  VSDRY: S-wave velocity of dry porous rock
%	  K0 MU0 RO0:  Mineral Bulk modulus shear modulus and density
%	  ROFLKFLNU: Pore fluid density bulk modulus and viscosity
%	  POR PERM: Porosity and absolute permeability of porous rock
%	  A: pore size parameter; usually about 1/6-1/7 of grain diameter
%	  ALFA: Tortuosity parameter （always > 1） usually between 1-3
%	  D1D2: frequency range from 10^D1 to 10^D2
%	  ‘PLOTOPT‘: optional input for plotting. If nothing is given plots of
%	  all six outputs versus frequency are displayed. The following options
%	  selects only one of the outputs for plotting.
%	  ‘PLOTOPT‘=‘VP1‘ displays VP1 （fast P-wave） dispersion curve
%		   =‘VP2‘ displays VP2 （slow P-wave） dispersion curve
%		   =‘VS‘  displays VS （shear wave） dispersion curve
%		   =‘QP1‘ displays P1 attenuation curve
%		   =‘QP2‘ displays P2 attenuation curve
%		   =‘QS‘  displays S attenuation curve
% Outputs: VP1: Fast P-wave velocities at all frequencies
%	   FREQ: Frequencies at which velocity and attenuation are calculated
%	   VP2: Slow P-wave velocities at all frequencies
%	   VS:  S-wave velocities
%	   Q1INV Q2INV QSINV: Fast and slow P-wave and S-wave attenuations
%
%	   See also BIOTHF BIOTHFB

% Written by T. Mukerji 1996

if nargin<15 % nargin  是matlab的内置函数，返回值为函数的输入参数的个数；
    opt=‘all6‘; 
end; 
ro=（1-por）*ro0+por*rofl;%平均法计算饱和岩石的密度；
rodry=（1-por）*ro0; %计算干岩心的密度；
omc=por*nu/（rofl*perm）;% 计算临界角频率；
mudry=rodry.*vsdry.^2; % 利用Vs_dry计算干岩心（骨架）的剪切模量；
kdry=rodry.*vpdry.^2-（4/3）*mudry; %利用Vp_dry计算干岩心（骨架）的体积模量；
d=k0*（1+por*（k0/kfl-1））; % 计算参数D，表达式见Handbook中P140；
h=kdry+（4/3）*mudry+（k0-kdry）.^2./（d-kdry）; % 计算参数H;
c=k0*（k0-kdry）./（d-kdry）; %计算参数C；
m=k0.^2./（d-kdry）; %计算参数M；
%从第47行到第58行都是为了计算q;
freq=logspace（d1d2100）; % 在[0^d1，10^d2]之间等对数间距的产生100个数据点，返回值为行向量，1*100；
om=2*pi*freq; %计算角频率om为行向量，长度为100;
zeta=sqrt（（rofl*a^2/nu）*om）;%计算中间参数Zeta，为行向量，1*100;
t=zeros（size（zeta））; %设置零矩阵，1*100；
idx=[zeta <= 1e3];% idx与Zeta具有相同的长度，如果Zeta（j）<=1e3则idx（j）=1;反之，idx（k）=0;
idxl=[zeta <= 1e-1];
t（idx）=exp（i*3*pi/4）*bessel（1exp（-i*pi/4）*zeta（idx））./...
        bessel（0exp（-i*pi/4）*zeta（idx））;
t（~idx）=（（1+i）/sqrt（2））*ones（sum（~idx）1）;
f=0.25*zeta.*t./（1+i*2*t./zeta）; % 计算参数F（Zeta）
f（idxl）=ones（sum（idxl）1）;
qf=alfa*rofl/por - （i*nu/perm）*f./om;% 计算参数q;
%%%%%%%%%%%%test of qf
%qf=alfa*rofl/por + （i*nu/perm）*f./om; this equation is wrong; sign is -
%%%%%%%%%%%

b0=rofl^2-ro*qf;
b1=h*qf+m*ro-2*c*rofl;
b2=c^2-m*h; %b0b1b2为求解Vp时用的参数；
sl1sqr=（-b1+sqrt（b1.^2-4*b2*b0））./（2*b2）;
sl2sqr=（-b1-sqrt（b1.^2-4*b2*b0））./（2*b2）;
slssqr=（ro*qf-rofl^2）./（mudry*qf）;
vp1=1./real（sqrt（s

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2013-07-02 15:43  biot 双向介质正演模拟matlab程序\
     文件     2221462  2013-04-23 18:01  biot 双向介质正演模拟matlab程序\250ms_snapshotz_f.dat
     文件     2220020  2013-04-23 18:01  biot 双向介质正演模拟matlab程序\250ms_snapshotz_s.dat
     文件        4104  2011-09-29 17:15  biot 双向介质正演模拟matlab程序\biot.m
     文件       10278  2011-12-22 16:01  biot 双向介质正演模拟matlab程序\biot_2_4_PML.m
     文件    13860613  2013-04-23 18:01  biot 双向介质正演模拟matlab程序\offset_vx_s
     文件    13860613  2013-04-23 18:02  biot 双向介质正演模拟matlab程序\offset_vz_s