这是一个matlab代码，用于执行Backus平均，以推导测井数据的低频模型。这一理论发表在Kumar（2013）上.zip

%------------------------------------------------------- % Matlab code to do Backus averaging as in Kumar (2013), % "Applying Backus averaging for deriving seismic anisotropy % of long-wavelength equivalent medium from well-log data" % Backus averaging of isotropic/VTI mixed layers % This program is used to compute data and plot Figure 3 % Author: Dhananjay Kumar, Email: dhananjaykumar@gmail.com % Input: layer properties in terms of Thomsen parameters (VTI) % Output: Thomsen's parameters (VTI) of the Effective Medium % How to run: save this file and type file.m on the matlab prompt % Workflow: % 1. Input parameters (need user input) % 2. Computation as described in paper % 3. Matlab plotting %------------------------------------------------------- clear all; %------------------------------------------------------- %1. Input parameters (Table 1) %------------------------------------------------------- no_rock = 2; %no of rock types present % rock type 1 - shale (VTI) vp1 = 3060; %m/s vs1 = 1490; %m/s rho1 = 2.42; %g/cc eps1 = 0.256; %Thomsen's parameters gam1 = 0.481; del1 = -0.051; % rock type 2 - sand (Isotropic) vp2 = 2950; %m/s vs2 = 1480; %m/s rho2 = 2.0; %g/cc eps2 = 0.0; %Thomsen's parameters gam2 = 0.0; del2 = 0.0; %------------------------------------------------------- % 2. Computation %------------------------------------------------------- %Layer's properties (Cij for VTI model) c33_1 = vp1*vp1*rho1; % rock type 1 c44_1 = vs1*vs1*rho1; c11_1 = (1+2*eps1)*c33_1; c66_1 = (1+2*gam1)*c44_1; c13_1 = sqrt(2*del1*c33_1*(c33_1-c44_1)+(c33_1-c44_1)^2)-c44_1; c33_2 = vp2*vp2*rho2; % rock type 2 c44_2 = vs2*vs2*rho2; c11_2 = (1+2*eps2)*c33_2; c66_2 = (1+2*gam2)*c44_2; c13_2 = sqrt(2*del2*c33_2*(c33_2-c44_2)+(c33_2-c44_2)^2)-c44_2; % Build one set of Cij matrix for all rock types for i = 1:no_rock if rem(i,2) == 1 %layer1 (shale) c11(i) = c11_1; c33(i) = c33_1; c44(i) = c44_1; c13(i) = c13_1; c66(i) = c66_1; elseif rem(i,2) == 0 %layer2 (sand) c11(i) = c11_2; c33(i) = c33_2; c44(i) = c44_2; c13(i) = c13_2; c66(i) = c66_2; end end j = 0; for rock1_fraction = 0:0.01:1 j = j+1; rock1_vol_fraction(j) = rock1_fraction; v(1) = rock1_fraction; %volume fraction: rock1 v(2) = 1-v(1); %volume fraction: rock2 c33_e = 1/sum((v./c33)); %Cij of effective model (eqns 8a-8e) c44_e = 1/sum((v./c44)); c13_e = sum(v.*(c13./c33))*c33_e; c66_e = sum(v.*c66); c11_e = sum(v.*c11)+c13_e*sum(v.*(c13./c33))-sum(v.*(c13.*c13./c33)); epsilon(j) = 0.5*(c11_e-c33_e)/c33_e; delta(j) = 0.5/(c33_e*(c33_e-c44_e))*((c13_e+c44_e)^2-(c33_e-c44_e)^2); gamma(j) = 0.5*(c66_e - c44_e)/c44_e; rho_e = v(1)*rho1 + v(2)*rho2; % Note that vertical Vp and Vs are independent of intrinsic VTI vp_e(j) = sqrt(c33_e/rho_e); %m/s vs_e(j) = sqrt(c44_e/rho_e); %m/s end % end of for loop %------------------------------------------------------- % 3. Plot output: Thomsen parameters as a function of rock1 volume %------------------------------------------------------- plot(rock1_vol_fraction,epsilon,'k.','linewidth',2); hold on; plot(rock1_vol_fraction,delta,'k-','linewidth',2); hold on; plot(rock1_vol_fraction,gamma,'k--','linewidth',2); hold on; ylim ([-0.1 0.5]); xlabel ('Rock type 1 (volume fraction)'); ylabel ('VTI Anisotropy'); legend ('epsilon','delta','gamma'); %-------------------------END---------------------------