zouxiang_i.m

沉陷变形计算。本程序可以给予概率积分模型计算沉陷变形的水平移动。

clear;
clc;
 syms x s z qy z2 xx yy tl zx tt;
 syms q d t qx_jf;
%  w=[0 0 0 0 0 0 0 0 0 0 0 0 0 0
%        0 -4 2 0 -5 7 2 1 7 -1 0 -10 1 -10
%        0 0 -4 1 8 -4 6 -3 9 7 8 -17 3 11
%        0 1 4 6 9 -2 9 -2 8 4 6 1 2 52
%        0 4 4 5 7 -3 9 -1 6 1 1 3 -9 46
%        0 15 22 28 35 18 27 11 12 10 11 4 -3 50	
%        0 16 13 19 26 13 20 6 9 9 9 2 -1 54	
%        0 21 25 35 32 13 16 4 7 7 9 3 1 58	
%        0 33 44 59 58 28 29 2 11 10 8 -4 4 62
%        0 40 62 93 95 53 39 7 15 14 7 1 -6 57
%        0 52 69 133 141 97 61 20 19 18 13 1 -1 55	
%        0 49 81 150 191 149 95 39 22 21 20 2 -2 52
%        0 46 65 158 223 228 174 75 45 40 30 10 0 47
%        0 38 59 139 217 243 197 95 55 55 26 16 -1 51
%        0 32 40 111 192 245 241 141 76 92 47 47 5 45
%        0 23 27 78 155 218 255 167 113 117 67 26 22 44
%        0 21 20 54 117 177 249 188 135 157 92 52 16 46
%        0 12 14 41 85 119 214 197 149 185 140 76 41 43			
%        0 9 6 25 61 72 156 175 168 204 185 101 66 39
%        0 7 7 31 36 44 111 143 160 215 222 150 82 44		
%        0 6 3 9 5 36 73 103 124 229 234 188 98 49];
%修正1506点下沉
   w=[0 0 0 0 0 0 0 0 0 0 0 0 0 0
       0 -4 2 0 -5 7 2 1 7 -1 0 -10 1 -10
       0 0 -4 1 8 -4 6 -3 9 7 8 -17 3 11
       0 1 4 6 9 -2 9 -2 8 4 6 1 2 52
       0 4 4 5 7 -3 9 -1 6 1 1 3 -9 46
       0 10 9 12 18 5 15 3 8 5 5 2 -4 50	
       0 16 13 19 26 13 20 6 9 9 9 2 -1 54	
       0 21 25 35 32 13 16 4 7 7 9 3 1 58	
       0 33 44 59 58 28 29 2 11 10 8 -4 4 62
       0 40 62 93 95 53 39 7 15 14 7 1 -6 57
       0 52 69 133 141 97 61 20 19 18 13 1 -1 55	
       0 49 81 150 191 149 95 39 22 21 20 2 -2 52
       0 46 65 158 223 228 174 75 45 40 30 10 0 47
       0 38 59 139 217 243 197 95 55 55 26 16 -1 51
       0 32 40 111 192 245 241 141 76 92 47 47 5 45
       0 23 27 78 155 218 255 167 113 117 67 26 22 44
       0 21 20 54 117 177 249 188 135 157 92 52 16 46
       0 12 14 41 85 119 214 197 149 185 140 76 41 43			
       0 9 6 25 61 72 156 175 168 204 185 101 66 39
       0 7 7 31 36 44 111 143 160 215 222 150 82 44		
       0 6 3 9 5 36 73 103 124 229 234 188 98 49];
[m,n]=size(w);
%计算累计下沉
   for i=1:1:m
    wl(i,1)=w(i,1);
    for j=2:1:n
        wl(i,j)=wl(i,j-1)+w(i,j);
    end
end
% 累计天数,以开采（12月26日起算）
d2=[0,10,21,30,42,61,81,102,118,136,156,177,204,241,501];
d3=[10,21,30,42,61,81,102,118,136,156,177,204,241,501];
%两时刻之间均值
d_j=[5,15.5,25.5,36,46.5,71,91.5,110,127,151,171.5,190.5,222.5,371];
%走向各点间距
s1=[0 93.8 140.3 172.9 201.7 236.4 266.2 295.6 329.2 358.3 390.3 420.3 474.8 504.3 533.6 566.6 593.7 624.3 658.2 685.4 715];
%以开切眼为原点。
s2=s1-390.3;
 %计算各点间实测倾斜
   for i=2:1:m
       for j=1:1:n
           w_w(i-1,j)=wl(i,j)-wl(i-1,j);
           qx(i-1,j)=w_w(i-1,j)/(s2(i)-s2(i-1));
       end
   end
  %两点点间平均位置
    for i=2:1:m
       s_pj0(i-1)=0; s_pj(i-1)=(s2(i)+s2(i-1))/2;
    end
     %三点之间平均距离
   for i=3:1:m
       s_k(i-2)=s2(i-2)+(s2(i)-s2(i-2))/2;
   end
    %计算实测曲率
   for i=3:1:m
       for j=1:1:n
           k(i-2,j)=(qx(i-1,j)-qx(i-2,j))/(s2(i)-s2(i-2));
       end
   end
  % 计算实测下沉速度
for i=1:1:m
    for j=2:1:15
       vs(i,j-1)=w(i,j-1)/(d2(j)-d2(j-1));
    end
end
% %倾斜图形
%   plot(s_pj,qx(:,1),'ro-',s_pj,qx(:,2),'b+-',s_pj,qx(:,3),'y*-',s_pj,qx(:,4),'k*-',s_pj,qx(:,5),'mx-')
%   xlabel('两点点间平均位置','fontweight','bold');ylabel('走向倾斜','fontweight','bold');
%   title('走向倾斜分布曲线1','fontsize',12,'fontweight','bold','fontname','隶书');figure
%   plot(s_pj,qx(:,6),s_pj,qx(:,7),'ro-',s_pj,qx(:,8),'b+-',s_pj,qx(:,9),'y*-',s_pj,qx(:,10),'k*-')
%     xlabel('两点点间平均位置','fontweight','bold');ylabel('走向倾斜','fontweight','bold');
%   title('走向倾斜分布曲线2','fontsize',12,'fontweight','bold','fontname','隶书');figure
%   plot(s_pj,qx(:,11),'cx-',s_pj,qx(:,12),'k*-',s_pj,qx(:,13),'ro-',s_pj,qx(:,14),'b+-')
%     xlabel('两点点间平均位置','fontweight','bold');ylabel('走向倾斜','fontweight','bold');
%   title('走向倾斜分布曲线3','fontsize',12,'fontweight','bold','fontname','隶书');figure
% %显示实测速度,距离
% plot(s2,vs(:,1),'k*-',s2,vs(:,2),'y*-',s2,vs(:,3),'b+-',s2,vs(:,4),'ro-')
% figure
% plot(s2,vs(:,5),'k*-',s2,vs(:,6),'y*-',s2,vs(:,7),'b+-',s2,vs(:,8),'ro-',s2,vs(:,9),'mx-')
% figure
% plot(s2,vs(:,10),'k*-',s2,vs(:,11),'y*-',s2,vs(:,12),'b+-',s2,vs(:,13),'ro-',s2,vs(:,14),'cx-')
%   xlabel('走向测点距离','fontweight','bold');ylabel('走向实测速度','fontweight','bold');
%title('走向下沉速度曲线','fontsize',12,'fontweight','bold','fontname','隶书');
% figure
% %显示实测速度,时间
% plot(d_j,vs(1,:),'k*-',d_j,vs(2,:),'y*-',d_j,vs(3,:),'b+-',d_j,vs(4,:),'ro-')
% figure
% plot(d_j,vs(5,:),'k*-',d_j,vs(6,:),'y*-',d_j,vs(7,:),'b+-',d_j,vs(8,:),'ro-')
% figure
% plot(d_j,vs(9,:),'k*-',d_j,vs(10,:),'y*-',d_j,vs(11,:),'b+-',d_j,vs(12,:),'ro-')
% figure
% plot(d_j,vs(12,:),'k*-',d_j,vs(13,:),'y*-',d_j,vs(14,:),'b+-',d_j,vs(15,:),'ro-')
% figure
% plot(d_j,vs(16,:),'b+-',d_j,vs(17,:),'ro-')
% figure
% plot(d_j,vs(18,:),'k*-',d_j,vs(19,:),'y*-',d_j,vs(20,:),'b+-',d_j,vs(21,:),'ro-')
%   xlabel('走向两次观测时间平均','fontweight','bold');ylabel('走向实测速度','fontweight','bold');
%   title('走向下沉速度曲线','fontsize',12,'fontweight','bold','fontname','隶书');figure
%实验
% yj=[0,50,100,130,150,200,250,300,500,1000];qys=[-200,-100,0,50,100,130,150,200,250,400];% xj=[0,50,100,130,150,200,250,500,800,1500];zxs=[-200,-100,0,100,200,300,500,800,1000,1200];
% yj=[0,189,378,567,756,2000];qys=[-189,0,94.5,189,383.5];
% for i=1:1:6
%     for j=1:1:5
%         y_jfs(i,j)=int(exp(-pi*(qys(j)-yy)^2/r^2),yy,0,yj(i));x_jfs(i,j)=int(exp(-pi*(zxs(j)-xx)^2/r^2),xx,0,xj(i));wgs(i,j)=w0/r^2*y_jfs(i,j)*x_jfs(i,j);
%     end
% end
% y_jfs=double(y_jfs);x_jfs=double(x_jfs);wgs=double(wgs);
%基于概率积分的各点下沉、下沉速度,w0采用采高1928，而不是最大下沉1340
%两时刻之间均值
w0=1928;r=189;
qy=60;zx=s2;xw=s2;lw=14;rq=163;
d_j=[5,15.5,25.5,36,46.5,71,91.5,110,127,151,171.5,190.5,222.5,371];
d3=[10,21,30,42,61,81,102,118,136,156,177,204,241,501];
kc=[11.2 44.6 70 104.5 125 148 196 246 286 316 371 433 515 606 606];
v=[1.9 2.2 2.8 3.1 2 2.5 2.5 2.5 2.6 2.8 3 3 2.5 2.5];
y_jf=w0/rq*int(exp(-pi*(qy-yy)^2/rq^2),yy,0,130);y_jf=double(y_jf);
for i=1:1:21
     for j=1:1:14
        %速度对结果影响很大
        wg(i,j)=y_jf/r*int(exp(-pi*(zx(i)-xx)^2/r^2),xx,0,kc(j));
        vg(i,j)=v(j)*y_jf/r*exp(-pi*(zx(i)-v(j)*d_j(j))^2/r^2);           
    end
end
vg=double(vg);wg=double(wg)
%概率倾斜
   for i=2:1:m
       for j=1:1:n
           w_wg(i-1,j)=wg(i,j)-wg(i-1,j);
           qxg(i-1,j)=w_wg(i-1,j)/(s2(i)-s2(i-1));
       end
   end
%取整
qxg=double(qxg);
%概率理论倾斜
for i=1:1:m-1
    for j=1:1:14        
    gqx(i,j)=y_jf/r*int(exp(-pi*(s_pj(i)-xx)^2/r^2)*(-2*pi*(s_pj(i)-xx)/r^2),xx,0,kc(j));
   end
end
gqx=double(gqx)
%由倾斜计算概率曲率
     for i=3:1:m
       for j=1:1:n
           kg(i-2,j)=(qxg(i-1,j)-qxg(i-2,j))/(s2(i)-s2(i-2));
       end
     end
%计算概率理论曲率
 for i=3:1:m
       for j=1:1:n
          gk(i-2,j)=(-2*pi/r^2)*y_jf/r*int(exp(-pi*(s_k(i-2)-xx)^2/r^2)*((-2*pi*(s_k(i-2)-xx)^2/r^2)+1),xx,0,kc(j));
       end
 end 
   gk=double(gk); 
% %显示理论速度，以距离为横轴
% plot(s2,vg(:,1),'ro-',s2,vg(:,2),'b+-',s2,vg(:,3),'y*-',s2,vg(:,4),'k*-',s2,vg(:,5),'mx-')
% figure
% plot(s2,vg(:,6),s2,vg(:,7),'ro-',s2,vg(:,8),'b+-',s2,vg(:,9),'y*-',s2,vg(:,10),'k*-')
% figure
% plot(s2,vg(:,11),'cx-',s2,vg(:,12),'k*-',s2,vg(:,13),'ro-',s2,vg(:,14),'k*-')
%   xlabel('走向测点距离','fontweight','bold');ylabel('走向概率速度','fontweight','bold');
%   title('走向下沉速度曲线','fontsize',12,'fontweight','bold','fontname','隶书');figure
% %显示理论速度，以时间为横轴
% plot(d_j,vg(1,:),'ro-',d_j,vg(2,:),'b+-',d_j,vg(3,:),'y*-',d_j,vg(4,:),'k*-')
% figure
% plot(d_j,vg(5,:),'mx-',d_j,vg(6,:),'ro-',d_j,vg(7,:),'b+-',d_j,vg(8,:),'y*-')
% figure
% plot(d_j,vg(9,:),'k*-',d_j,vg(10,:),'cx-',d_j,vg(11,:),'k*-',d_j,vg(12,:),'ro-',d_j,vg(13,:),'b+-',d_j,vg(14,:),'y*-')
% figure
% plot(d_j,vg(15,:),'k*-',d_j,vg(16,:),'cx-',d_j,vg(17,:),'k*-',d_j,vg(18,:),'ro-',d_j,vg(19,:),'cx-',d_j,vg(20,:),'k*-',d_j,vg(21,:),'ro-')
%   xlabel('走向两次观测时间平均','fontweight','bold');ylabel('走向概率速度','fontweight','bold');
%   title('走向下沉速度曲线','fontsize',12,'fontweight','bold','fontname','隶书');figure
%吴侃曲线  
 syms xx yy;
%以开切眼为原点。
xw=s2;yw=60;rw=189;lw=14;hw=400;gw=v*2/rw;
y_jfk=w0/rq*int(exp(-pi*(yw-yy+lw)^2/rq^2),yy,0,130);
%li=400*ctg(88);
for i=1:1:21
    for j=1:1:14 
       fy=v(j)*d_j(j);
    ww(i,j)=y_jfk/r*int(exp(-pi*(xw(i)-xx)^2/rw^2)*(1-exp(-gw(j)*(d3(j)-xx/v(j)))),xx,0,kc(j));
    vw(i,j)=gw(j)*y_jfk/r*int(exp(-pi*(xw(i)-xx)^2/rw^2)*exp(-gw(j)*(d_j(j)-xx/v(j))),xx,0,fy);
   end
end
ww=double(ww);vw=double(vw)
%吴侃理论倾斜1，直接理论计算
for i=1:1:m-1
    for j=1:1:14 
    kqx(i,j)=y_jfk/r*int(exp(-pi*(s_pj(i)-xx)^2/rw^2)*(1-exp(-gw(j)*(d3(j)-xx/v(j))))*(-2*pi*(s_pj(i)-xx)/rw^2),xx,0,kc(j));
   end
end
kqx=double(kqx)
%吴侃理论倾斜2,以下沉为基础计算倾斜
   for i=2:1:m
       for j=1:1:n
           w_wk(i-1,j)=ww(i,j)-ww(i-1,j);
           qxk(i-1,j)=w_wk(i-1,j)/(s2(i)-s2(i-1));
       end
   end
%取整
qxk=double(qxk)
%吴侃理论曲率
     for i=1:1:m-2
       for j=1:1:n
           kk(i,j)=(-2*pi/rw^2)*y_jfk/r*int(exp(-pi*(s_k(i)-xx)^2/rw^2)*(1-exp(-gw(j)*(d3(j)-xx/v(j))))*(1-2*pi*(s_k(i)-xx)^2/rw^2),xx,0,kc(j));
       end
   end
kk=double(kk);
%倾斜实测、概率、吴侃比较，以距离为横轴
for i=1:1:n
    plot(s_pj,gqx(:,i),'ro-',s_pj,qx(:,i),'b+-',s_pj,kqx(:,i),'y*-')
    text(100,-30,num2str(i));figure
end
% for i=1:1:n
%     plot(s_pj,qxg(:,i),'ro-',s_pj,gqx(:,i),'b+-')
%     text(100,-30,num2str(i));figure
% end
% for i=1:1:n
%     plot(s_pj,qxk(:,i),'ro-',s_pj,kqx(:,i),'b+-')
%     text(100,-30,num2str(i));figure
% end

% %由吴侃倾斜计算的曲率
     for i=3:1:m
       for j=1:1:n
           kkw(i-2,j)=(qxk(i-1,j)-qxk(i-2,j))/(s2(i)-s2(i-2));
       end
   end
% %以时间为横轴
% for i=1:1:m-1
%     plot(d_j,qxg(i,:),'ro-',d_j,qx(i,:),'b+-',d_j,qxk(i,:),'y*-',d_j,l_qx(i,:),'k*-')
%     text(200,-30,num2str(i));figure
% end

% %曲率实测、概率、吴侃比较，以距离为横轴
% for i=1:1:n
%     plot(s_k,gk(:,i),'ro-',s_k,k(:,i),'b+-',s_k,kk(:,i),'y*-',s_k,l_k(:,i),'k*-',s_k,kg(:,i),'cx-')
%     text(100,-30,num2str(i));figure
% end
% % 曲率实测、概率、吴侃比较，以时间为横轴
% for i=1:1:m-2
%     plot(d_j,gk(i,:),'ro-',d_j,k(i,:),'b+-',d_j,kk(i,:),'y*-',d_j,l_k(i,:),'k*-')
%     text(100,-30,num2str(i));figure
% end
%%显示吴侃速度,距离
% plot(s2,vw(:,1),'ro-',s2,vw(:,2),'b+-',s2,vw(:,3),'y*-',s2,vw(:,4),'k*-')
% figure
% plot(s2,vw(:,5),'mx-',s2,vw(:,6),'ro-',s2,vw(:,7),'b+-',s2,vw(:,8),'y*-',s2,vw(:,9),'k*-')
% figure
% plot(s2,vw(:,10),'cx-',s2,vw(:,11),'k*-',s2,vw(:,12),'ro-',s2,vw(:,13),'b+-',s2,vw(:,14),'y*-')
%   xlabel('走向测点距离','fontweight','bold');ylabel('走向吴侃速度','fontweight','bold');
%   title('走向下沉速度曲线','fontsize',12,'fontweight','bold','fontname','隶书');figure
%显示吴侃速度,时间
% plot(d_j,vw(1,:),'ro-',d_j,vw(2,:),'b+-',d_j,vw(3,:),'y*-',d_j,vw(4,:),'k*-')
% figure
% plot(d_j,vw(5,:),'mx-',d_j,vw(6,:),'ro-',d_j,vw(7,:),'b+-',d_j,vw(8,:),'y*-',d_j,vw(9,:),'k*-',d_j,vw(10,:),'cx-')
% figure
% plot(d_j,vw(11,:),'k*-',d_j,vw(12,:),'ro-',d_j,vw(13,:),'b+-',d_j,vw(14,:),'y*-',d_j,vw(15,:),'k*-',d_j,vw(16,:),'cx-',d_j,vw(17,:),'k*-')
% figure
% plot(d_j,vw(18,:),'ro-',d_j,vw(19,:),'b+-',d_j,vw(20,:),'y*-',d_j,vw(21,:),'ro-')
%   xlabel('走向两次观测时间平均','fontweight','bold');ylabel('走向吴侃速度','fontweight','bold');
%   title('走向下沉速度曲线','fontsize',12,'fontweight','bold','fontname','隶书');
% %速度实测、概率、吴侃比较，以时间为横轴
% for i=1:1:21
%     plot(d_j,vg(i,:),'ro-',d_j,vs(i,:),'b+-',d_j,vw(i,:),'y*-')
%     text(100,0,num2str(i));figure
% end
% %速度以距离为横轴
% for i=1:1:14
%     plot(s2,vg(:,i),'ro-',s2,vs(:,i),'b+-',s2,vw(:,i),'y*-',s2,zeros(size(s2)),'k')
%     text(100,0,num2str(i));figure
% end
% %下沉实测、概率、吴侃比较，以时间为横轴
% for i=1:1:21
%     plot(d3,-wg(i,:),'ro-',d3,-wl(i,:),'b+-',d3,-ww(i,:),'y*-')
%     text(100,-30,num2str(i));figure
% end
% %以距离为横轴
% for i=1:1:14
%     plot(s2,-wg(:,i),'ro-',s2,-wl(:,i),'b+-',s2,-ww(:,i),'y*-')
%     text(200,-30,num2str(i))
%     figure
% end
% 克诺特曲线
% g=0.2;
% g2=1-exp(-g*tt);
% d=1:1:100;
% ezplot(g2,[0 100]);
%   xlabel('时间','fontweight','bold');ylabel('速度','fontweight','bold');
%   title('克诺特曲线','fontsize',12,'fontweight','bold','fontname','隶书');figure
%for i=1:1:21
    %ezplot(g1(i),[0 200]);
%end
%王金庄速度曲线
% vmax=13;l=45;a=16;
% for i =1:1:21
%     v(i)=vmax/(1+((s1(i)+l)/a)^2);
% end
% plot(s1,v,'ro-')
% figure
% v=vmax/(1+((xx+l)/a)^2);
% ezplot(v,[-100 40]);
% figure
%计算实测和概率下沉差
for i=1:1:14
    qx_gl(i)=0;qx_kl(i)=0;
end
  for i=1:1:m-1
      for j=1:1:14
          qx_gxd(i,j)=int32((qx(i,j)-gqx(i,j))/(qx(i,j)+0.0001)*100);
          qx_g(i,j)=qx(i,j)-gqx(i,j);
          qx_g_pf(i,j)=qx_g(i,j)^2; 
                    
          qx_k(i,j)=qx(i,j)-kqx(i,j);
          qx_kxd(i,j)=int32((qx(i,j)-kqx(i,j))/(qx(i,j)+0.01)*100); 
          qx_k_pf(i,j)=qx_k(i,j)^2;
      end
      qx_g_h(i)=sqrt(sum(qx_g_pf(i,:))/13);
      qx_k_h(i)=sqrt(sum(qx_k_pf(i,:))/13);
  end
  %w_gxd=double(w_gxd);w_kxd=double(w_kxd);
  for i=1:1:14
        qx_g_l(i)=sqrt(sum(qx_g_pf(:,i))/13);
        qx_k_l(i)=sqrt(sum(qx_k_pf(:,i))/13);
  end
  %出差前进行到比较各种参数下概率、吴侃预计与实测下沉的误差比较，回来后可考虑调整r、qy和w0等参数。