unction main()
%%%%%%%%%%%%%%%%%%%%%%%-------注入与产出数据处理-------%%%%%%%%%%%%%%%%%%%%%%
close all; clear; clc;
%形成生产井数据
[num,txt]=xlsread('pro.xls');
proname = num(:,1); %n*1---井名的列向量
num = num(:,2:end); %num(n,:)为第n口井在各时间点上的数据,行为井
%num = num';
[hNum,lNum] = size(num);
for i=1:2:lNum
pro1(:,(1+i)/2) = num(:,i) + num(:,i+1); %产液
pro1_hsl(:,(1+i)/2) = num(:,i+1)./(num(:,i) + num(:,i+1)+0.00001*ones(length(num(:,i+1)),1)); %含水率
end
%%%%形成注入数据
[num,txt]=xlsread('inj.xls');
injname = num(2:end,1)'; %n*1---井名的行向量
num = num';
inj1=num(:,2:end); %inj1(:,n)为第n口井在各时间点上的数据,行为时间
%xlswrite('coord_pro.xls',proname'); %xlswrite('coord_inj.xls',injname')%形成头文件
% xlswrite('连通系数.xls', proname', '连通系数', 'A2') ;
% xlswrite('连通系数.xls', injname, '连通系数', 'B1') ;
% xlswrite('连通系数.xls', xlsnada, '连通系数', 'B2') ;
%提取计算数据,3个月为1组
pro1 = pro1'; pro1_hsl = pro1_hsl';
for i=1:3:size(pro1,1)
pro((2+i)/3,:) = pro1(i,:)+pro1(i+1,:)+pro1(i+2,:);
pro_hsl(:,(2+i)/3) = (pro1_hsl(i,:)+pro1_hsl(i+1,:)+pro1_hsl(i+2,:))./3; %含水率
inj((2+i)/3,:) = inj1(i,:)+inj1(i+1,:)+inj1(i+2,:);
end
%clear pro1_hsl; clear pro1;
%xlswrite('inj3.xls',inj); xlswrite('pro3.xls',pro) %形成文件
%%%%%%%%%%%%%%%%%%%%%%-------井位坐标处理-------%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[cdpro,tmp] = xlsread('coord_pro.xls'); %cdpro(n,:)--第n口井xy坐标
procdname=cdpro(:,1); %procdname列井名
[cdinj,tmp] = xlsread('coord_inj.xls');%cdinj(n,:)--第n口井xy坐标,
injcdname = cdinj(:,1); %列井名
%计算注入井和各生产井之间的距离
for i = 1:size(cdinj,1)
for j=1:size(cdpro,1)
inj_pro_dis(j,i) = norm(cdpro(j,:)-cdinj(i,:));
end
end
for i=1:size(inj_pro_dis,2)
[y I]=sort(inj_pro_dis(:,i),'ascend');
dis_sort(:,i)=I;
end
% for i=1:size(dis_sort,1)
% aa = reshape(dis_sort(1:i,:),1,size(dis_sort,2)*i);
% size(unique(aa))
% end
% clear aa;
%%%%%%%%%%%%%%%%%%%%%-------灰色关联度分析-------%%%%%%%%%%%%%%%%%%%%%%%%%%%
%xi=pro_hsl; [rows,cols]=size(xi); p= 0.3 %含水率波动太小
xi=pro'; [rows,cols]=size(xi); p= 0.3
for k = 1:size(inj,2)
x0 = inj(:,k)'; %母数据
%求差值
for i = 1:rows
deta0i(i,:)=abs(xi(i,:) - x0);
end
detamax = max(max(deta0i)); detamin = min(min(deta0i)); %得到最大最小值
%求灰色关联系数
for i = 1:rows %rows个油井
for j = 1:cols
kc0i(i,j) = (detamin + p*detamax)/(deta0i(i,j)+p*detamax);
end
end
out(:,k) = mean(kc0i,2);
%井距离限制
for i = 1:rows
if inj_pro_dis(i,k) > inj_pro_dis(dis_sort(12,k),k)
out(i,k) = out(i,k) - 1;
end
% if(size(find(xi(i,:)<=0.01),2) > size(xi,2)*0.7) %去掉很多数据为0的井
% out(i,k) = out(i,k) - 2;
% end
end
end
for i=1:size(out,2)
[y I]=sort(out(:,i),'descend');
well(:,i)=I;
end
% for i=1:size(well,1)
% aa = reshape(well(1:i,:),1,size(well,2)*i);
% size(unique(aa))
% end
% clear aa;
clear inj_pro_dis; clear dis_sort;clear out; clear xi; clear pro_hsl;
% pro=pro1;
% inj=inj1;
%%%%%%%%%%%%%%%%%%%%%%----拟合----%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nadaxls = zeros(size(proname,2),size(injname,2));
for i = 1:size(inj,2) %注入井
%for i = 1:2 %注入井
pron = [pro(:,well(1,i)) pro(:,well(2,i)) pro(:,well(3,i)) pro(:,well(4,i)) pro(:,well(5,i)) ...,
pro(:,well(6,i)) pro(:,well(7,i)) pro(:,well(8,i)) pro(:,well(9,i)) pro(:,well(10,i))];
thiswell = [proname(well(1,i)) proname(well(2,i)) proname(well(3,i)) proname(well(4,i)) proname(well(5,i)) ...,
proname(well(6,i)) proname(well(7,i)) proname(well(8,i)) proname(well(9,i)) proname(well(10,i))];
j = 1;
injn = inj(:,i);
while(1)
if sum(pron(j,:) < 0.01) == size(pron,2)
pron(j ,:) = [];
injn(j ,:) = [];
j = j-1;
end
j = j+1;
if(j == size(pron,1))
break
end
end
[nadaout(:,i),dout(:,i)] = cal_nada(pron,injn,i,thiswell);
i
for kk = 1:10
if (nadaout(kk,i) < 0.01)
nadaout(kk,i) = 0;
end
nadaxls(well(kk,i),i) = nadaout(kk,i);
delayxls(well(kk,i),i) = dout(kk,i);
end
end
dlmwrite('xlsnada',nadaxls);
dlmwrite('xlsdelay',delayxls);
dlmwrite('xlsnada8',nadaout);
dlmwrite('xlsdelay8',dout);
dlmwrite('nadaout',nadaout);
dlmwrite('nadout',dout);
dlmwrite('bbb',nadaxls);
sum(sum(nadaxls')>0)
nadaout = dlmread('xlsnada8');
dout = dlmread('xlsdelay8');
nadaxls = dlmread('xlsnada');
%%%%%%%%%%%%%%%%%%%-----绘图-----%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:size(well,2)
figure(i);
plot(inj1(:,i),'bv-','LineWidth',2);hold on;
plot(pro1(:,well(1,i)),'r-','LineWidth',1);hold on;
plot(pro1(:,well(2,i)),'y-','LineWidth',1);hold on;
plot(pro1(:,well(3,i)),'g-','LineWidth',1);hold on;
plot(pro1(:,well(4,i)),'c-','LineWidth',1);hold on;
plot(pro1(:,well(5,i)),'k-','LineWidth',1);
s1 = num2str(injname(i)); s2 = num2str(proname(well(1,i)));
s3 = num2str(proname(well(2,i))); s4 = num2str(proname(well(3,i)));
s5 = num2str(proname(well(4,i))); s6 = num2str(proname(well(5,i)));
legend(s1,s2,s3,s4,s5,s6);
Title(strcat('注入井:',num2str(injname(i))));
s1 = char('z_');
s1 = strcat(s1,num2str(injname(i)));
s1 =strcat(s1,'.png');
% s1 =strcat(char(injname(i)),'.png');
print(gcf,'-dpng',s1)
close(figure(gcf)) ;
end
%成果图
%画井位,写井位名称
figure(size(well,2) + 1);
plot(cdpro(:,1),cdpro(:,2),'ro','MarkerFaceColor','r');hold on;
text(cdpro(:,1)+5,cdpro(:,2)+5, num2str(procdname)); hold on;
plot(cdinj(:,1),cdinj(:,2),'bo','MarkerFaceColor','b');hold on;
text(cdinj(:,1)+5,cdinj(:,2)+5, num2str(injcdname)); hold on;
mindis = 500;
for i=1:size(well,2)
plot_nada(cdinj(i,:),cdpro(well(1,i),:),nadaout(1,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(2,i),:),nadaout(2,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(3,i),:),nadaout(3,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(4,i),:),nadaout(4,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(5,i),:),nadaout(5,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(6,i),:),nadaout(6,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(7,i),:),nadaout(7,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(8,i),:),nadaout(8,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(9,i),:),nadaout(9,i),mindis);hold on;
plot_nada(cdinj(i,:),cdpro(well(10,i),:),nadaout(10,i),mindis);hold on;
end
axis square;
print(gcf,'-dpng','成果图.png');
nadaout
end
%注入点坐标,生产井坐标,连通系数(0-1),三角形高度(最小距离的一半)
function plot_nada(injpos,propos,nada,heigh)
if nada > 0.001
distmp = norm(injpos-propos);
fdk = nada*heigh/(distmp - nada*heigh);
B = (injpos + fdk*propos )/(1+fdk);
l = nada * heigh * tan(nada/2);
a = injpos(1); b = injpos(2);
m = B(1); n = B(2);
C1x = m - l*(a-m)/sqrt((a - m)^2 + (b - n)^2);
C2x = m + l*(a-m)/sqrt((a - m)^2 + (b - n)^2);;
C1y = n - (b - n)*(C1x - m)/(a-m);
C2y = n - (b - n)*(C2x - m)/(a-m);
fill([a C1x C2x a],[b C1y C2y b],'r');
end
end
详细详细每一步都解释我真的看不懂
最新发布