数模-城市表层土壤重金属污染分析

题目

数据

问题一的求解

用MATLAB软件对所给数据进行插值拟合得出调查区的地形图和8种主要重金属元素在该城区的空间分布图，再用MATLAB软件对所给数据进行分析得出功能区散点图：

Matlab griddata函数功能介绍

城区地形分布图 代码：

A=xlsread('cumcm2011A附件_数据.xls',1,'A4:E322');
x=A(:,2);y=A(:,3);z=A(:,4);
scatter(x,y,5,z)%散点图 5是散点的大小
title('功能区')
figure
[X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值
pcolor(X,Y,Z)
shading interp %伪彩色图
title('功能区')
figure,contourf(X,Y,Z) %等高线图
title('功能区')
figure,contour(X,Y,Z) 
title('功能区')
figure
surf(X,Y,Z)%三维曲面
title('功能区')

功能区分布散点图 代码：

A=xlsread('cumcm2011A附件_数据.xls',1,'A4:E322');
x=A(:,2);y=A(:,3);z=A(:,5);
x1=find(z==1);
x_1=x(x1(:));
y_1=y(x1(:));
scatter(x_1,y_1,20,'d');
hold on
x2=find(z==2);
x_2=x(x2(:));
y_2=y(x2(:));
scatter(x_2,y_2,20,'h');
hold on
x3=find(z==3);
x_3=x(x3(:));
y_3=y(x3(:));
scatter(x_3,y_3,20,'s');
hold on
x4=find(z==4);
x_4=x(x4(:));
y_4=y(x4(:));
scatter(x_4,y_4,20,'p');
hold on
x5=find(z==5);
x_5=x(x5(:));
y_5=y(x5(:));
scatter(x_5,y_5,20,'x');
title('功能区分布')
legend('生活区','工业区','山区','主干道区','公园绿地区')

重金属在该城区空间分布图的代码：

A=xlsread('cumcm2011A附件_数据.xls',1,'A4:E322');
B=xlsread('cumcm2011A附件_数据.xls',2,'B4:I322');
x=A(:,2);y=A(:,3);
for k=1:8
    z=B(:,k);   
    [X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值
    figure
    surf(X,Y,Z)%三维曲面
end

不同的重金属的空间分布图（Z轴是浓度大小）：

论文中的效果图：

前几问可以看相关论文

城市表层土壤重金属污染分析

第三问

代码：

%砷、镉、铬、铜...的污染源分布
A=xlsread('cumcm2011A附件_数据.xls',1,'A4:E322');
B=xlsread('cumcm2011A附件_数据.xls',2,'B4:I322');
ss={'As','Cd','Cr','Cu','Hg','Ni','Pd','Zn'};
x=A(:,2);y=A(:,3);%X轴 Y轴
for k=1:8
    z=B(:,k);
    [X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值
    z=Z';
    for i=2:99
        for j=2:99
            if (z(i,j)>z(i-1,j))&&(z(i,j)>z(i+1,j))&&(z(i,j)>z(i,j+1))&&...
                    (z(i,j)>z(i,j-1))&&(z(i,j)>z(i-1,j-1))&&(z(i,j)>z(i-1,j+1))...
                    &&(z(i,j)>z(i+1,j-1))&&(z(i,j)>z(i+1,j+1));
                z(i,j)=1000;
            end;
        end;
    end;
    [ii,jj]=find(z==1000);
    disp(ii');disp(jj');
    subplot(2,4,k)
    scatter(ii,jj,'*')
    title(ss{k})
end

筛选污染源的代码（以砷为例，其它元素的程序只需调整参数即可）

A=xlsread('cumcm2011A附件_数据.xls',1,'A4:E322');
B=xlsread('cumcm2011A附件_数据.xls',2,'B4:I322');
ss={'As','Cd','Cr','Cu','Hg','Ni','Pd','Zn'};
x=A(:,2);y=A(:,3);%X轴 Y轴
z=B(:,1);
[X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值
z=Z';
k=z;
for i=2:99
    for j=2:99
        if (z(i,j)>z(i-1,j))&&(z(i,j)>z(i+1,j))&&(z(i,j)>z(i,j+1))&&...
                (z(i,j)>z(i,j-1))&&(z(i,j)>z(i-1,j-1))&&(z(i,j)>z(i-1,j+1))...
                &&(z(i,j)>z(i+1,j-1))&&(z(i,j)>z(i+1,j+1));
            z(i,j)=1000;
        end;
    end;
end;
[ii,jj]=find(z==1000);
c=[ii,jj];
n=size(ii,1);
for i=1:n
    kk(i)=k(ii(i),jj(i));%这里的kk就是对应的插值出来的浓度Z的搜索值
end
pp=kk';
kkk=find(kk>15);%As的一级污染源 返回的是索引
kkk=kkk';
n=size(kkk,1);
for i=1:n
    rr(i)=pp(kkk(i));
end
rr %返回的是一级污染源的浓度
for i=1:n
    cc(i,1)=c(kkk(i),1);%c是坐标
    cc(i,2)=c(kkk(i),2);
end
cc %返回的是一级污染源的坐标
scatter(cc(:,1),cc(:,2),'*')
axis([0,100,0,100])
title(ss{1})

rr返回的是一级污染源的浓度
cc返回的是一级污染源的坐标

转载请注明出处，欢迎讨论和交流! 【文章转自：美国站群服务器 http://www.558idc.com/mgzq.html处的文章,转载请说明出处】