2015年认证杯SPSSPRO杯数学建模
C题 荒漠区动植物关系的研究
原题再现：
环境与发展是当今世界所普遍关注的重大问题, 随着全球与区域经济的迅猛发展, 人类也正以前所未有的规模和强度影响着环境、改变着环境, 使全球的生命支持系统受到了严重创伤, 出现了全球变暖、生物多样性消失、环境污染等全球性的环境问题, 并已经严重影响到了全球人类社会的发展。在探讨环境与发展问题的过程中, 人们越来越认识到了现代社会发展过程中自然——社会——经济复合生态系统的复杂性, 以及生态学理论在解决这些问题中的重要性。
干旱区是全球生态系统中的重要类型之一, 也是目前全球开发较晚的区域之一, 因此, 积极开展干旱区的生态学理论与实践研究, 对于干旱区当前面临的重大环境问题的解决, 以及未来防患于未然的科学决策均具有极其重要的现实意义。作为我国三大自然区域之一的西北干旱区, 由于其大规模、高强度的开发历史较短, 因此, 与其它区域相比较而言, 其境内蕴藏了丰富的待开发自然资源, 也奠定了其在我国未来经济建设中的举足轻重的战略地位, 并担负着重要的历史使命, 西部大开发战略的实施即是最显著的证明。因此, 积极开展和深化干旱区的生态学研究, 对于该区域的经济发展与生态环境保护具有深远的理论意义与实践价值。
生态研究与资源利用是分不开的, 荒漠区是我国典型的温带荒漠和干旱脆弱生态系统, 生态环境条件十分严酷, 动物的可利用资源在数量和质量上与湿润区、半干旱区存在差异, 啮齿动物的分布具有明显的区域性特征。由于近年来人为干扰不断加重, 使得该地区的荒漠化日益严重。依赖于植物生存的动物种群和群落格局随之受到了明显影响。
啮齿动物群落是荒漠生态系统食物链上必不可少的消费者, 对荒漠的利用与保护有至关重要作用。许多物种群体与人的干扰具有密切关系, 干扰的一个突出作用是导致生态系统中各类资源的改变和生态系统结构的重组, 导致异质性环境的形成。有关不同干扰方式下, 栖息地破碎化过程中研究群落的变化特征是当前景观生态学和群落生态学研究的前沿。
第二阶段问题：
1. 请结合附件一和附件二的数据，建立合理的数学模型，来评估由人类活动造成的荒漠地区生态退化的程度。
2. 分析一个荒漠地区处于半退化、退化等不同阶段时，是否可以通过减少人为干扰，或者采用补充人工植被的方法来促使该地区的生态环境恢复正常。如果可行，请给出量化的实施方案；如果不可行，请指出造成这种不可逆性的原因。
整体求解过程概述(摘要)
目前，土地的荒漠化已经成为一个世界热点话题。土地的荒漠化会导致土地生产力的下降、土地资源的丧失以及对环境的破坏，对我们人类的发展与生存产成极大的影响。因此，研究对在荒漠地域生态退化程度的判断标准以及在此标准下在各阶段的补救措施具有重要的现实意义。
对于问题一，首先，根据分析可知，生态退化的程度可以在调查的动植物的种群状态上反映出来，故我们将植物因子与啮齿动物优势种百夹捕获率作为生态退化程度的评价指标，将附件中的重复项作为样本。同时我们发现数据中存在着部分的缺失，由于数据的时序关系较弱，插值等数据补充方法难以应用，故将存在缺失数据的样本删除。接着，为了减少数据偶然性的影响，我们利用基于共享型最近邻居相似度算法（SNN）进行了异常值的筛检，共去除了9个异常样本（见表 2）。然后，对数据进行归一化后，我们建立了投影寻踪模型，利用基于实数编码的加速遗传算法求解得到最佳投影方向为a*=(0.0761 ,0.1167 ,0.1735 ,0.1113 ,0.2068 ,0.3520 ,0.2208 ,0.1362 , 0.5483 ,0.5081 ,0.3811)。利用最佳投影方向计算各样本的投影值（见附录3），由于投影值能很好地表现高维数据的分类结构特征，故将其作为评价值进行评价模型的构建。最后，为获得生态退化、半退化与未退化状态的区分标准，我们采用了K-Means聚类算法将样本的评价值聚为3类，分析评价值的属性，以生态状态来解释各分类类别，并以类间的中间距离作为不同状态的分界线得到区分标准。据此讨论不同人类活动造成的荒漠地区生态退化情况（见表 8），发现开垦会造成土地极大的破坏，在开垦的影响下，88.37%的地域处于生态退化状态，过牧对土地的负面影响次之，轮牧对土地的负面影响较于前两者而言相对较弱一些。
对于问题二，为了直观地反映过牧、轮牧、开垦三种荒漠地区的半退化和退化的动态变化情况，我们采用元胞自动机（CA）进行仿真模拟。在不改变当前人为干扰和人工植被的情况下，三种荒漠区呈现出不同的生态退化情况，轮牧区较轻，过牧区略重，开垦区非常严重。然后，分别考虑减少人为干扰和增加人工植被两种方法进行治理，结果发现轮牧区和过牧区在减少35%人为干扰或增加一倍人工植被的情况下呈现出良好的恢复趋势，而开垦区减少65%人为干扰也难有起色，只有增加3倍人工植被才能慢慢恢复。说明对于过牧区和轮牧区能有效治理，而开垦区的治理花费太大，需要寻找新的方法。最后我们针对模型做了敏感性分析，并进行了中肯的评价和适当的推广。
问题分析：
对于问题一，首先分析附件中的数据是否存在异常值或缺失值，并进行数据处理。对于缺失的数据，由于附件中数据的时序关系等较弱，故我们将其去除不予考虑。因为数据的规模较小，存在偶然性的可能较大，为去除偶然性对评估的影响，考虑到样本数据为高维形式，我们采用基于共享型最近邻居相似度算法（SNN）进行异常值的筛检。然后，对数据进行归一化后，我们建立了投影寻踪模型，由于投影值可以很好地表现高维数据的分类结构特征，我们以此来获取样本的投影值作为评价值，最后采用 K-Means算法聚类进行样本的评估。
对于问题二， 我们选用元胞自动机（CA）分别对轮牧区、过牧区、开垦区进行仿真模拟，观察三个地区在未来几年内非退化、半退化、退化三种区域面积的动态变化情况，然后根据题目要求，修改元胞演化参数，分别模拟减少人为干扰和增加人工植被两种治理措施，观察此时三个地区的变化情况，并对最终结果的非退化、半退化、退化面积进行统计分析，以便能得出较好的量化治理措施。
模型假设：
1) 假设附件中同一月的重复项为在西北地区同一个月在不同小地域调查所得；
2) 由于评估生态退化的指标多种多样，指标的选取影响着评估模型的效果，假设题中所给的几项可作为生态退化的典型指标；
3) 假设剔除异常值之后的数据都具有统计意义，且数据来源真实可靠；
4) 假设附件中的小区域数据能有效地反映当地整体荒漠化的情况。
论文缩略图：
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
dis=distance(Vari_F);
[s,d]=rk(dis);
snn=rs(20,s,d);
[c,p]=cho(15,snn);
[n,m]=size(Vari_F);
Variable=zeros(n-p,m);
p1=1;
for i=1:n
k=1;% 指示变量
for j=1:c
if i==c(1,j)
k=0;
end;
end;
if k==0
continue;
else
Variable(p1,:)=Vari_F(i,:);
p1=p1+1;
end;
end;
%%
function a=distance(x2)
[n,m]=size(x2);
a=zeros(n,n);
c=0;
for i=1:n
for j=i:n
for k=1:m
b(k)=x2(i,k)-x2(j,k);
end;
for k=1:m
c=c+b(k)^2;
end;
a(i,j)=sqrt(c);
c=0;
end;
end;
for i=1:n
for j=1:n
a(j,i)=a(i,j);
end;
end;
%%
function [a,d]=rk(dis)% a为序号矩阵，d为距离矩阵
[n,m]=size(dis);
a=zeros(n,m);
for i=1:n
for j=1:m
a(i,j)=j;
end;
end;
for i=1:n
for k=1:m-1
for j=1:m-k
if dis(i,j)>dis(i,j+1)
t=a(i,j);
a(i,j)=a(i,j+1);
a(i,j+1)=t;
t=dis(i,j);
dis(i,j)=dis(i,j+1);
dis(i,j+1)=t;
end;
end;
end;
end;
d=dis;
%%
function [c,p]=cho(t,snn)% j为相似度矩阵,c为编号，p为总个数
[n,m]=size(snn);
for i=1:n
for j=1:m
if i==j
snn(i,j)=0;
end;
end;
end;
c=zeros(1,n);
for i=1:n
k(i)=1;
end;
for i=1:n
for j=1:m
if snn(i,j)>t
k(i)=0;
k(j)=0;
end;
end;
end;
p=0;
for i=1:n
if k(i)~=0
c(1,p+1)=i;
p=p+1;
end;
end;
d=[];e=[];
X=Variable;
[m,n]=size(X);
for k=1:50
x=normalize(X);
N=400;M=10;Ci=7;n=11;DaiNo=2;ads=1;
[a1,b1,ee,ff]=RAGA(x,N,n,Pc,Pm,M,DaiNo,Ci,ads);
d=[d,a1];e=[b1;e];
end
[a2 b2]=max(d);
e1=e(b2,:);
ff=e1*x';
%%
function a=normalize(b)
D=b;
maxi=max(D);
mini=min(D);
[n m]=size(D);
for i=1:n
for j=1:m
D(i,j)=(D(i,j)-mini(j))/(maxi(j)-mini(j));
end
end
a=D;
%%
function [a,b,mmin,mmax]=RAGA(xx,N,n,M,DaiNo,Ci,ads)
if ads==0
ad='ascend';
else
ad='descend';
end
mm1=zeros(1,n);mm2=ones(1,n);
Pc=0.8;Pm=0.2;
for z=1:Ci
while 1==1
for p=1:n
bb(p)=unifrnd(mm1(p),mm2(p));
end
temp=sum(bb.^2);
a=sqrt(bb.^2/temp);
y=Feasibility(a);
if y==1
v(i,:)=a;
break;
end
end
end
for s=1:DaiNo
for i=1:N
fv(i)=Target(xx,v(i,:));
end
[fv,i]=sort(fv,ad);
v=v(i,:);
arfa=0.05;
q(1)=0;
for i=2:N+1
q(i)=q(i-1)+arfa*(1-arfa)^(i-2);
end
for i=1:N
r=unifrnd(0,q(N+1));
for j=2:N+1
if r>q(j-1) & r<=q(j)
vtemp1(i,:)=v(j-1,:);
end
end
end
while 1==1
CrossNo=0;
v1=vtemp1;
for i=1:N
r1=unifrnd(0,1);
if r1 < Pc
CrossNo=CrossNo+1;
vtemp2(CrossNo,:)=v1(i,:);
v1(i,:)=zeros(1,n);
end
end
if CrossNo~=0 & mod(CrossNo,2)==0
break;
elseif CrossNo==0|mod(CrossNo,2)==1
vtemp2=[];
end
end
shengyuNo=0;
for i=1:N
if v1(i,:)~=zeros(1,n)
shengyuNo=shengyuNo+1;
vtemp3(shengyuNo,:)=v1(i,:);
end
end
for i=1:CrossNo
r2=ceil(unifrnd(0,1)*(CrossNo-i+1));
vtemp4(i,:)=vtemp2(r2,:);
vtemp2(r2,:)=[];
end
for i=1:2:(CrossNo-1)
while 1==1
r3=unifrnd(0,1);
v20(i,:)=r3*vtemp4(i,:)+(1-r3)*vtemp4(i+1,:);
v20(i+1,:)=(1-r3)*vtemp4(i,:)+r3*vtemp4(i+1,:);
temp1=sum(v20(i,:).^2);
temp2=sum(v20(i+1,:).^2);
v2(i,:)=sqrt(v20(i,:).^2/temp1);
v2(i+1,:)=sqrt(v20(i+1,:).^2/temp2);
if Feasibility(v2(i,:))==1 & Feasibility(v2(i+1,:))==1
break;
end
end
end
v3=[vtemp3;v2];
while 1==1
MutationNo=0;
v4=v3;
for i=1:N
r4=unifrnd(0,1);
if r4<Pm
MutationNo=MutationNo+1;
vtemp5(MutationNo,:)=v4(i,:);
v4(i,:)=zeros(1,n);
end
end
if MutationNo~=0
break;
end
end
shengyuNo1=0;
for i=1:N
if v4(i,:)~=zeros(1,n)
shengyuNo1=shengyuNo1+1;
vtemp6(shengyuNo1,:)=v4(i,:);
end
end
DirectionV=unifrnd(-1,1,1,n);
for i=1:MutationNo
tempNo=0;
while 1==1
tempNo=tempNo+1;
v5(i,:)=sqrt(((vtemp5(i,:)+M*DirectionV).^2)./sum((vtemp5(i,:)
+M*DirectionV).^2));
y=Feasibility(v5(i,:));
if tempNo==200
v5(i,:)=vtemp5(i,:);
break;
elseif y==1
break;
end
M=unifrnd(0,M);
end
end
vk=[v5;vtemp6];
v=vk;
end
for i=1:N;
fv(i)=Target(xx,v(i,:));
end
[fv,i]=sort(fv,ad);
v=v(i,:);
vk=v;
vv=vk(1:20,:);
t=1:n;
mm1(t)=min(vv(:,t));
mm2(t)=max(vv(:,t));
mmin(z,:)=mm1;
mmax(z,:)=mm2;
if abs(mm1-mm2)<=0.00001
break;
end
end
a=fv(1);
b=vv(1,:);
%%
function y=Feasibility(a)
b=sum(a.^2);
if abs(b-1)<=0.00001
y=1;
else
y=0;
end
%%
function y=Target(x,a)
[m,n]=size(x);
for i=1:m
s1=0;
for j=1:n
s1=s1+a(j)*x(i,j);
end
z(i)=s1;
end
Sz=std(z);
R=0.1*Sz;
s3=0;
for i=1:m
for j=1:m
r=abs(z(i)-z(j));
t=R-r;
if t>=0
u=1;
else
u=0;
end
s3=s3+t*u;
end
end
Dz=s3;
y=Sz*Dz;