2016年第五届数学建模国际赛小美赛
A题 臭氧消耗预测
原题再现：
臭氧消耗包括自1970年代后期以来观察到的若干现象：地球平流层（臭氧层）臭氧总量稳步下降，以及地球极地附近平流层臭氧（称为臭氧空洞）春季减少幅度更大。除了这些众所周知的平流层现象外，还有春季极地对流层臭氧消耗事件。
据认为，臭氧消耗的主要原因是含氯源气体的存在，包括氟氯化碳和相关卤代烃、氧化亚氮等。1985年，20个国家签署了《保护臭氧层维也纳公约》。1987年，43个国家的代表签署了《蒙特利尔议定书》。在蒙特利尔，与会者同意将氟氯化碳的生产冻结在1986年的水平，并到1999年将生产减少50%。
自从通过和加强《蒙特利尔议定书》导致氟氯化碳排放量减少以来，大气中最重要化合物的浓度一直在下降。请建立一个数学模型来预测未来50年的臭氧水平。
整体求解过程概述(摘要)
空气是生命的基础，空气质量决定着人们的生活质量。众所周知，平流层臭氧层的减少防止了紫外线对生命形式的损害。如何准确地预测臭氧浓度对我们具有重要意义。
首先，详细介绍了臭氧产生的原因、臭氧空洞和臭氧损耗。然后对臭氧的时空分布进行分析，为预报做准备。分析结果表明，臭氧总量随纬度的增加而增加，也随季节的变化而变化。具体地说，我们的研究表明，北极地区春季臭氧含量显著增加，冬季臭氧含量显著减少。南极洲正好相反。
其次，考虑到全球臭氧浓度具有季节性和周期性，建立时间序列模型对臭氧柱浓度序列数据进行拟合。收集了南极和北极的臭氧最小值数据，建立了南极和北极的ARMA模型。经过分析，建立了南极is-ARMA（5,7）模型，建立了北极is-ARMA（11,3）模型。根据模型拟合结果和方程，得到未来50年臭氧最小值的预测值。
第三，考虑到平流层臭氧消耗物质，如氯、含溴碳氢化合物和氮氧化物。它们对臭氧在平流层中的分解具有催化作用，对平流层臭氧的稳定性构成威胁。为此，分析了大气中某些消耗臭氧物质对臭氧密度的影响，建立了灰色多元回归模型，剔除了不必要的ODS变量。建立多元回归方程，最终预测未来臭氧消耗物质的水平，并在此基础上预测未来50年空气中臭氧浓度。
本文还分析了模型的误差，分析了模型的优缺点，并提出了改进和完善的建议。特别是对模型的不足之处进行了扩展和改进。
模型假设：
我们正在寻找的数据是准确的或有微小的误差，但不影响整个模型。
人们的行为意识在未来50年内，变化规律与1979-2015年相似或模型没有显著影响。
不管紫外线对臭氧浓度的影响。
问题分析：
国内外科学家对大气臭氧变化趋势做出了重要贡献。Donald研究了北半球O3总量的长期变化，结果表明，O3总量在20世纪60年代初达到最低值，然后开始上升。20世纪70年代初O3总量达到最大值，50年代以来一直呈下降趋势。利用1975～1990年的TOMS数据，Stolarski利用统计模型计算了各纬度地区臭氧年平均趋势，发现随着纬度的变化，赤道附近的趋势几乎为零，逐渐向两极增加，臭氧减少主要集中在冬春季。
魏定文分析了北京和昆明1979-1993年地面观测资料，发现两地臭氧总量不断减少，尤其是1991-1993年。周秀吉根据卫星资料发现了青藏高原臭氧“槽”区。YanLiu探索了一种预测臭氧含量演变的分析方法，将臭氧总量划分为3～5个状态，用状态间的转移频率代替状态间的转移概率，最大频率状态为臭氧即将向趋势方向演化，并举例说明，效果较好。
但这一问题尚未完全解决，而现有文献或多或少存在一些不完善之处，需要我们加以完善。
我们的工作
将臭氧浓度分为臭氧柱浓度、臭氧总量和臭氧密度。
简要介绍了臭氧产生的原因、臭氧空洞和臭氧损耗。
首先，我们分析了臭氧的时空分布，为下一步的预测做准备。
其次，由于春季平流层臭氧减少，春季也有极地对流层臭氧消耗事件。我们假设全球臭氧浓度具有季节性和周期性，因此我们考虑一个时间序列模型来拟合臭氧柱浓度序列。
由于臭氧层主要发生在南极和北极地区，收集了南极和北极的臭氧最小值数据，分别设计了南极和北极的ARMA模型，并对模型进行了拟合和预测。
认为人类活动直接或间接地向平流层提供一些消耗臭氧层物质（ODS），如氯、含溴碳氢化合物和氮氧化物等，对平流层臭氧的分解具有催化作用，导致平流层臭氧的稳定性受到威胁。
为此，分析了大气中主要消耗臭氧层物质对臭氧层密度的影响，建立了灰色多元回归模型，并在此基础上对未来50年的臭氧层密度进行了预测。
最后对模型进行了误差分析，并提出了模型的优缺点。提出了改进和完善模型的建议。针对模型的不足，对模型进行了扩展和改进。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
clc,clear
a=[230 379 242 271 140 69 13 0.084120283
234 382 240 266 150 71 14 0.086837264
237 382 238 250 161 72 16 0.086521226
239 381 236 225 171 73 18 0.090599057
242 380 234 200 177 73 20 0.085811321
245 378 232 172 184 73 22 0.08509434
246 376 230 145 190 73 24 0.090320755
248 374 228 121 196 72 26 0.088084906
249 372 225 101 201 72 28 0.082063679
249 369 222 84 204 71 30 0.083731132
250 366 220 70 207 71 32 0.087933962
250 363 217 58 208 70 33 0.082316038
250 360 215 48 211 69 35 0.075070755
249 356 213 40 214 69 36 0.080023585
249 353 210 34 214 68 37 0.079304245
248 350 207 28 214 68 39 0.07925
247 347 205 24 212 67 40 0.075245283
246 344 201 20 210 66 42 0.069929245
244 341 199 17 208 66 43 0.073382075
243 338 196 14 206 65 45 0.074410377
242 335 193 12 203 65 46 0.075929245
241 332 191 10 201 64 48 0.067299528
239 330 189 8 197 63 49 0.068033019
238 329 186 7 193 63 50 0.068995283];
x1=a(:,1);x2=a(:,2);x3=a(:,3);x4=a(:,4);
x5=a(:,5);x6=a(:,6);x7=a(:,7);y=a(:,8);n=size(a,1);
% To produce CFC-12; there was no significant correlation between halons and annual mean ozone density
subplot(3,3,1);plot(x1,y,'*'); title('CFC-12')
subplot(3,3,2);plot(x2,y,'*');title('CFC-11')
subplot(3,3,3);plot(x3,y,'*');title('CCl4')
subplot(3,3,4);plot(x4,y,'*');title('CH3CCl3')
subplot(3,3,5);plot(x5,y,'*');title('halons')
subplot(3,3,6);plot(x6,y,'*');title('CFC-113')
subplot(3,3,7);plot(x7,y,'*');title('HCFCs')
% For multiple regression
x=[ones(n,1) x2 x3 x4 x6 x7];
[b,bint,r,rint,stats]=regress(y,x);
b,bint,stats,
%k=5,n=24,F(5,18)=4.96<29.7546 , established
%p<0.0001，established
%Gray prediction GM (1,1), prediction of CFC-11, CCl4, CH3CCl3, CFC-113, HCFCs for the next 50
years 0
% For the next 50 years for the annual mean ozone density
A=[326.7551524 185.5355533 6.380777381 63.26204622 57.79727923
324.2989558 183.3726964 5.38926077 62.82179442 60.73986646
321.8612222 181.2350527 4.551817109 62.38460641 63.8322673
319.4418129 179.1223284 3.844504817 61.95046087 67.0821091
317.0405902 177.0342328 3.247102627 61.51933663 70.49740753
314.6574173 174.9704789 2.742531475 61.09121266 74.08658636
312.2921585 172.930783 2.316366237 60.66606807 77.8584982
309.9446792 170.9148646 1.956423323 60.24388215 81.82244642
307.6148458 168.9224465 1.652412369 59.82463429 85.988208
305.3025256 166.9532548 1.395642039 59.40830405 90.36605771
303.0075869 165.0070186 1.178771557 58.99487113 94.96679343
300.7298991 163.0834705 0.995600839 58.58431536 99.80176277
298.4693326 161.1823458 0.840893236 58.17661672 104.8828911
296.2257586 159.3033833 0.710225833 57.77175532 110.2227109
293.9990493 157.4463246 0.599862993 57.36971143 115.8343927
291.7890781 155.6109142 0.506649567 56.97046544 121.7317776
289.5957191 153.7969 0.427920686 56.57399787 127.9294114
287.4188475 152.0040323 0.36142558 56.18028938 134.4425805
285.2583392 150.2320648 0.305263228 55.78932079 141.2873494
283.1140714 148.4807538 0.257828011 55.40107301 148.4806006
280.9859218 146.7498585 0.217763809 55.01552712 156.0400762
278.8737694 145.0391408 0.183925231 54.63266431 163.9844213
276.777494 143.3483656 0.15534487 54.25246591 172.3332306
274.696976 141.6773005 0.131205645 53.87491338 181.1070963
272.6320972 140.0257155 0.110817443 53.4999883 190.3276588
270.58274 138.3933837 0.093597388 53.1276724 200.0176605
268.5487876 136.7800806 0.079053179 52.7579475 210.2010017
266.5301243 135.1855843 0.066769011 52.39079559 220.9027993
264.5266352 133.6096757 0.056393695 52.02619875 232.149449
262.5382062 132.0521381 0.047630611 51.66413921 243.9686905
260.564724 130.5127572 0.040229234 51.3045993 256.3896757
258.6060764 128.9913215 0.033977965 50.9475615 269.4430407
256.6621518 127.4876217 0.028698088 50.59300838 283.1609814
254.7328396 126.0014511 0.024238658 50.24092267 297.5773328
252.8180299 124.5326054 0.020472184 49.89128718 312.7276525
250.9176137 123.0808825 0.017290987 49.54408486 328.6493084
249.0314827 121.6460829 0.01460412 49.19929879 345.3815711
247.1595297 120.2280093 0.012334768 48.85691214 362.9657102
245.3016481 118.8264667 0.010418054 48.51690823 381.4450967
243.457732 117.4412624 0.00879918 48.17927046 400.8653096
241.6276765 116.072206 0.007431865 47.84398237 421.2742484
239.8113774 114.7191091 0.006277019 47.51102762 442.7222515
238.0087313 113.3817858 0.005301625 47.18038995 465.2622197
236.2196356 112.0600522 0.004477799 46.85205325 488.9497476
234.4439884 110.7537265 0.003781989 46.5260015 513.8432599
232.6816886 109.4626292 0.003194301 46.20221881 540.0041559
230.9326359 108.1865826 0.002697935 45.88068937 567.4969609
229.1967307 106.9254114 0.002278699 45.56139752 596.3894854
227.4738742 105.6789421 0.00192461 45.24432767 626.7529922
225.7639683 104.4470034 0.001625542 44.92946437 658.6623722];
B=[ones(50,1) A];
C=B*b;
%CFC-11 Gray prediction
clc,clear
x0=[379 382 382 381 380 378 376 374 372 369 366 363 360 356 353 350 347 344 341 338
335 332 330 329];% Enter the raw data
x0=x0.^(1/2);
n=length(x0);
A(1,1)=exp(-2/(n+1));A(1,2)=exp(2/(n+1));
lamda=x0(1:n-1)./x0(2:n); % Calculate the level ratio
range=minmax(lamda);
A(2,:)=range;
x1=cumsum(x0);
for i=2:n
z(i)=0.5*(x1(i)+x1(i-1)); % Calculate the Immediate Mean Generation Number
end
B=[-z(2:n)',ones(n-1,1)];
Y=x0(2:n)';
u=B\Y;
x=dsolve('Dx+a*x=b','x(0)=x0');
x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});
yuce1=subs(x,'t',[0:n-1]);
yuce2=subs(x,'t',[0:n+49]);
digits(6),y=vpa(x);
yuce=[x0(1),diff(yuce1)];
yuce3=[x0(1),diff(yuce2)];
epsilon=x0-yuce ; % The residuals are calculated
delta=abs(epsilon./x0); % Calculate the relative error
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda;% Ratio deviation
yuce=yuce.^(2);
yuce3=yuce3.^(2);
