寻找污染源

%寻找污染源；clcclearA=load('at1.txt');B=load('at2.txt');b=B(:,2:9);c=[3.6 130 31 13.2 35 12.3 31 69];for i=1:length(b)p(i,:)=b(i,:)./c;e(i)=max(p(i,:));w=mean(p(i,:));d(i)=sqrt((e(i)^2+w^2)/2);endw1=1;r1=[];A1=[];p1=[];w2=1;r2=[];A2=[];p2=[];w3=1;r3=[];A3=[];p3=[];w4=1;r4=[];A4=[];p4=[];w5=1;r5=[];A5=[];p5=[];for i=1:length(b)switch A(i,5)case 1r1(w1)=d(i);A1(w1,:)=A(i,:);p1(w1,:)=p(i,:);w1=w1+1;case 2r2(w2)=d(i);A2(w2,:)=A(i,:);p2(w2,:)=p(i,:);w2=w2+1;case 3r3(w3)=d(i);A3(w3,:)=A(i,:);p3(w3,:)=p(i,:);w3=w3+1;case 4r4(w4)=d(i);A4(w4,:)=A(i,:);p4(w4,:)=p(i,:);w4=w4+1;otherwiser5(w5)=d(i);A5(w5,:)=A(i,:);p5(w5,:)=p(i,:);w5=w5+1;endendfor i=1:8As=find(p(:,1)>5);Cd=find(p(:,2)>5);Cr=find(p(:,3)>5);Cu=find(p(:,4)>10);Hg=find(p(:,5)>10);Ni=find(p(:,6)>5);Pb=find(p(:,7)>5);Zn=find(p(:,8)>5);endfigure x=A(:,2);y=A(:,3);z=A(:,4);[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x),200)',linspace(min(y),max(y),200),'v4');%插值C=contour(X,Y,Z,9);clabel(C) %二维图;axis square hold onplot(A(Hg,2),A(Hg,3),'bO')%浓度差较大的样本点ww1=1;rr1=[];AA1=[];pp1=[];ww2=1;rr2=[];AA2=[];pp2=[];ww3=1;rr3=[];AA3=[];pp3=[];ww4=1;rr4=[];AA4=[];pp4=[];ww5=1;rr5=[];AA5=[];pp5=[];for i=1:length(Hg)switch A(Hg(i),5)case 2rr2(ww2)=d(Hg(i));AA2(ww2,:)=A(Hg(i),:);pp2(ww2,:)=p(Hg(i),:);ww2=ww2+1;case 4rr4(ww4)=d(Hg(i));AA4(ww4,:)=A(Hg(i),:);pp4(ww4,:)=p(Hg(i),:);ww4=ww4+1;otherwiserr5(ww5)=d(Hg(i));AA5(ww5,:)=A(Hg(i),:);pp5(ww5,:)=p(Hg(i),:);ww5=ww5+1;endendhold onplot(AA2(:,2),AA2(:,3),'r*')hold onplot(AA4(:,2),AA4(:,3),'r+')w1=p(Hg(1:5),2)/sum(p(Hg(1:5),2));x1=w1'*A(Hg(1:5),2);y1=w1'*A(Hg(1:5),3);hold onplot(x1,y1,'bh','MarkerSize',16)text(x1,y1,'污染源');w2=p(Hg(14:16),2)/sum(p(Hg(14:16),2));x2=w2'*A(Hg(14:16),2);y2=w2'*A(Hg(14:16),3);hold onplot(x2,y2,'bh','MarkerSize',16)text(x2,y2,'污染源');w3=p(Hg(6:8),2)/sum(p(Hg(6:8),2));x3=w3'*A(Hg(6:8),2);y3=w3'*A(Hg(6:8),3);hold onplot(x3,y3,'bh','MarkerSize',16)text(x3,y3,'污染源');w4=p(Hg(11:12),2)/sum(p(Hg(11:12),2));x4=w4'*A(Hg(11:12),2);y4=w4'*A(Hg(11:12),3);hold onplot(x4,y4,'bh','MarkerSize',16)text(x4,y4,'污染源');title('Hg元素的污染源位置')[x1 y1;x2 y2;x3 y3;x4 y4]