2011高教社杯全国生数学建模竞赛一等奖题城市表层土壤重金属污染分析

先假设各重金属毒性及其它性质相同，运用公式

Pijj平均值的比值作为金属污染程度，再运用CiCij 7 三维地形 时空结 地质环境演 影响因子权

MATLAB ( 度

i区域中调查重金属的点 土样干密度（g/cm3 i区域中j重金属元素污染 度

XS

某区域各金属因子j在k

D(j）扩散系数，v为孔隙流体的

xf

粒子受到的浮 粒子受到的流体阻

将生活区、工业区、山区、交通区、公园绿地区用i表示。i=(1、2、3、4、5)Pij表示，Asμg/g)、Cd(ng/g)、Cr(μg/g)、(μg/g)、Hg(ng/g)Ni(μg/g)Pb(μg/g)Zn(μg/g)j=(1、2、…、

n

JI)下图（图1附件：AsCdCr9CuHg8NiPb6Zn设重金属元素污染程度Cij

jcij即 Ci

C2C3C4C5即 该城区内不同区域重金属的污染程NR式中：NRD

若对论域（研究的范围）U中的任一元素x，都有一个数A（x）∈[0，1]与之对应，则AU上的模糊集，A（x）称为xA的隶属度。当xU中变动时，A（x）就是一个函数，称为A的隶属函数。隶属度A（x）越接近于1，表示x属于A的程度越高，A（x）越接近于0表示x属于A的程度越低。用取值于区间[0，1]的隶属函数A（x）表征x属于A的程度高低例如A(x)=表示模糊集“年老”的隶属函数，A表示模糊集“年老”，当年龄x≤50A（x）=0表明x不属于模糊集A（即“年老），当x≥100时，A（x）=1表明x完全属于A，当50x〈1000〈A（x〈1x越接近100，A（x）1，x属于A的程度就越高。这样的表达方法显然比简单地说：“100100岁以下的人就不年老。”更为合理a、金属因子隶属度D线。设土壤环境质量分为m个级别,则 V=(1,2,⋯,m)

(PSjk jk

Sjk1S (SjkPSjk1(PSjkXjkSjk—某区域各金属因子j在k级指标(j1,2,⋯,n;k1,2,⋯m —各金属因子的实测浓度 D

b、金属因子权重向量R的确定:ffCj

nmn

j1式中,xjkjfjj j个重金属污染物的权重值,且Cj1rj R=[r1、r2,r3,⋯,rn ：本文所选用的八种元素，选用国家《土壤环境质量标准》[18](GB15618表（4） PH=6.5～7.5PHAsCdCrCuHgNiPbZn根据Hakanson：Zn=l，Cr=2Cu=5Pb=5N=5，As=10Cd=30Hg=40代表毒性越大。按照（7）、（8）两式，将表（1）的数据和毒性指数代入计算，得各区域各重金属权重值，其中Rr1、r2,r3⋯r8[0.09620.28380.06800.05853723、0.0211、0.0563、 N2N3= N4结合从整体的角度分析泥沙颗粒在水溶液中对重金属污染物的吸附遵循兰格茂尔吸附运动学方程dNKc(bNKNd dtNkc(bN)k (Bhu)By uuugy

u2

1

(Bhu)By(hs)B(ss

(Bhu)By(hs)By0 ABhA3By0及L=2hB 1(AE'c)1NA31NLsN1usN1E'sN1tuxA 1 1t t x 1xx NNNN'其中，N=N0，对（11） N1ksc(bsN)ks 2 m ' h=常数，u常数，00,sk( ), u E 1 1 1 N1kc(bN)k 2b=0.543g/kg,k1=0.0076k2=0.000841/sc

c0

x式中：c0 (x)= t0

x0,

0,

t0,

' tuxE1x20,ct0c0(x),cx0c0(t),c

于泥沙吸附需要时间,因此河段上泥沙吸附量的峰值总是滞后于水相浓度的峰值。散作用,峰型变化相同条件下重金属污染物输移,尽管图形相似,但重金属污染物的峰值因泥沙吸附而下降很多,特别是后半部减小更为明显,且随时间增加明显加大。Ft(j)F(j)e(j)/ )(

孔隙中的溶质的质量浓度 — —土样干密度（gcm3 vdFt(j)(AB)( AebssbB ．Ft(j)((j))n， R(j)(D((j))(j))v(j) 度j而变化的非线性的扩散系数，v(bn()1n/e)Ft(j) 1nnFt(j)) 1nnFt( x(DFt( Fj( 其中：t Ft0(j i1/ i1/ x RddFi(j) (j)Fi1(j)Fi(j) (j)Fi i1/ i1/ x Fi(j)Fi1(xi

— dwGF

(t 由式（24）得 ) w 其 分别为气体和 2r

子密度，通常《。令 g

w w(t)g(1et)etw(0)test0将w(t w(t)g(1et)

w(t)w(t)w(t)

t0

由上式求出w(t)的方差w2(t)，由于w(t)和 不相关，则

w2(t)w2(0)e2t)2K(1e2tw(0)w(t)w2(0)et

wwww

K22由（26）w(ttg(1e(tt)e(tt)w(0)(ttestt0w(0)etw(t)g(1et)test w(ttg(1e(tt)e(tt)w(t(ttestt0 将最后一项离散化得：w(tt)g(1t)tw(t) w(t)g(1et)etw(0) 方差为：w2(t)e2tw2(0) 22K

2t)2

2t

tw(t)

t)0 ) )u u(tt)uetu(t)u

t0 v(tt)vetv(t)

e

)05 w(tt)wg(1 ) w(t)w(1 )05有题可知，结合计算精度和时间的考虑，时间步长t(3—5)r DV InC(t1)其中，C(t)为t时刻的平均浓度，C(t)为t时刻的平均浓度Dt2

（图六c.在本文中，由于重金属的传播有固、液、气三种相属，而且，这三种相属又是相互联系，相互转换的，不可很明确的区分开来，故，我们将地壳运金属含生活工绿化程植被净 业污 度化同时还需要这个地区历年的地质环境影响因素资料及同一时期6yij表示第i城区地质环境的影响因素j的参数值。（假设，四个参数所对应的权重分别为ljj1,2,3,4)gi(i1,2,3,4,56)素的下属因素能直接查出，同上方法可以通过已知信息计算出yij。 lj(j1,2,3,4)=gi(i1,2,3,4,5,6)*可得各影响因素的权重值ljj1,2,3,4glj

进行权重值ljj1,2,3,4)二．缺点杨启航,数学建模,北京:高等教育出版社实用计算器4.50版软件, 2011年9月9 解可新，最优化方法，天津：天津大学出版社MATLAB %只需从excel输入三列数据，格式为：[xi %采样点坐标x值y=d(:,2);% 采样点坐标y值z=d(:,3);% shadinginterpholdonforfori=1:44fori=45:80fori=81:146fori=147:284fori=285:319%采样点坐标x%采样点坐标y%d=data;%只需从exc