backstrip.m

我自己编的计算地质埋藏史用的一维回拨程序。

clear all;
clc;
t=[0,55,65,100,140,160,210];                                      %   地层年代(Ma)
y0=[0,2116,2721,3689,4148,4148,5200];                            %   地层资料(m)
fai0=[0.63,0.49,0.70,0.63,0.63,0.56,0];                          %   地表孔隙度（今到古）                          %  地表孔隙度
c=[0.51,0.27,0.71,0.51,0.51,0.39,0];                             %   压实系数(1/km)
c=c*(1e-3);
lay=7;                                                              %   地层层数＋1
row=1e+3;                                                           %   水密度
rom=3.33e+3;                                                        %   地幔密度(kg/m3)
rosg=[0,2.72,2.65,2.71,2.72,2.71,2.68];                          %   各层骨架密度
rosg=rosg*(1e+3);
diet=[0,-80,-100,-150,-150,-150,-150];                           %   各时期的古海平面(从今到古）
wd=[172,220,150,0,0,0,0];                                          %   各时期的古水深                                                        %   地层深度为现今水深补偿
y0=y0-wd(1);
y0(1)=0;
D=1e+24;                                                             %   抗挠强度(Nm)
lamda=400*(1e+3);                                                   %   盆地宽度
yy=y0(lay)-y0;                                                       %   总沉降计算
for i=1:lay                                                          %   从古到今
    j=lay-i+1;
    [yyy(j),m]=depth(fai0(j),c(j),y0,lay,i);                      %   去压实 
    [y(j),ros,yyyy]=eqb(m,i,fai0(j),c(j),rosg,rom,lay);              %   沉积物负荷校正
    bendc=bendcoeff(rom,ros,D,lamda);                              %   此时挠曲系数求取
    y1(j)=bendc*(y(j)-diet(i)*row/(rom-row))+(wd(i)-diet(i));    %   真正的构造沉降                          
    m
end
figure(1);
plot(t,yy,'r',t,y,'g',t,yyy,'b',t,y1,'y');
legend('去压实前','负荷校正后','去压实后','真正的构造沉降');
xlabel('年代'),ylabel('地层厚度'),title(' 构造沉降图 ')
%--------------------压实深度与时间关系图------------------------------
figure(2);
for  k=1:lay  
   for i=k:lay
       j=lay-i+1;
       [yyy(j),m]=depth(fai0(j),c(j),y0,lay,i);
       y2(j)=m(i-k+1);
   end
   for i=1:lay-k+1
       tt(i)=t(i);
   end
   plot(tt,y2);
   hold on;
end
hold off;
%--------------------脱压实参数表---------------------------------------------