2020年华为杯数学建模竞赛E题论文和代码

大雾演化规律的量化分析与预测模型研究

随着人工智能技术的进步，近年来基于视频的路况能见度检测方法有了较大的创新，该方法可以克服传统激光能见度仪耗资巨大和探测范围小等问题，然而现有方法大多是只选取少量视频、截取图像中的某些固有特征进行能见度数值的间接估计。事实上，大雾的形成和消散有其自身的规律，视频资料中蕴含了丰富的大雾变化过程信息，充分利用这些信息，不仅可以提高能见度估计精度，也可以对大雾的消散进行预测。这对我国的高速公路管理和航空飞行具有重要的意义。

基于此，本文旨在通过收集和分析与能见度相关的气象因子信息和图像特征信息，评估各变量特征与能见度之间的相关性程度，并以此为基础在构建地面气象观测与能见度之间的多变量非线性回归模型；同时，充分利用视频的连续信息建立基于视频数据的能见度估计深度学习模型。然后利用图像中不同景深的物体在不同能见度下的亮度差异构建不依赖实测数据的能见度估计模型；最后运用时间序列分析模型量化能见度随时间的变化规律，建立数学模型预测大雾变化趋势以及其消散时间。

针对问题一，要求建立地面气象观测数据与能见度的具体关系式。为此，需要寻找出影响能见度的气象因子，从而建立能见度和气象因子的多元非线性回归模型。研究中结合文献调研，初步选取地面气象观测数据中可能影响能见度的因子，进而运用相关性分析评估各因子之间的内在关联，最终确定回归模型的变量因子。然后运用SPSS和Matlab软件对历史实测数据进行时空变化趋势做统计分析研究能见度与各气象因子的数学关系，以此构建出不同气候模式下估计能见度的多元非线性回归模型；在求解参数时，利用遗传算法得到多变量非线性回归模型的最优解并检验参数合理性。结果表明，利用该方法可以从定量和定性两个角度建立能见度和气象因子的回归方程，具有一定的可靠性和正确性。

针对问题二，要求建立基于视频数据与实测能见度数据的深度学习反演模型。为此，需要寻找出影响能见度的图像特征因子，从而构建能见度与视频数据的深度学习模型。首先，通过已有研究初步选取影响能见度的图像特征，然后运用相关性分析筛选最能影响雾天能见度的图像特征向量以及图像特征的窗格大小，同时评估各窗格区域内的特征值与能见度的关系，选取最优兴趣窗口，进而构建回归模型的参数向量；然后以多元线性回归方程建立参数向量与能见度之间的统计关系，并利用 BP 神经网络对线性拟合的残差进行估计修正多元线性回归的预测结果。结果表明，该模型能够同时兼顾线性和非线性的特点，反演精度较高，可以实现对能见度的精准估计。

针对问题三，要求建立非直接测量能见度的能见度估计算法。为此，需要获取图片深度信息估计能见度。从图像数据来看，近点有较清晰的道路信息，所以可借助道路线进行梯度检测求得边缘信息；已知车道线的长度 6m，利用亮度对比法结合光学模型的传播可以求出摄像机参数；远点因为大雾而看不到道路的尽头，因此我们利用求定的相机参数来获取离相机最远可分辨的像素，即人眼可分辨的亮度差异；然后使用暗原色先验原理获得目标点的透射率，求出大气折光系数，进而对能见度进行估计。结果表明其中 2016 年 4月14日6点30分至7点39分时段能见度变化趋势较缓，集中在47m附近，波动范围±10m。

针对问题四，要求根据能见度随时间的变化规律，建立数学模型预测大雾变化趋势以

及何时消散。为此，需要对能见度变化序列进行时序分析，主要包括倾向变化率、平稳性、

周期性、偏相关和自相关，进而通过线性回归模型、灰色预测模型和ARIMA模型模拟大雾的变化趋势以及消散过程。结果表明，问题三等到的能见度随时间变化具有一定的周期

性，且整体上保持恒定，上下波动不大；线性回归模型的结果显示下一时刻的大雾变化趋势逐渐减弱，且在8点50分56秒大雾开始消散即达到能见度为150m以上；ARIMA时序预测模型和灰色预测模型的结果显示大雾的变化趋势在短时间内逐渐加重，但速率较缓。

关键词：BP 神经网络、多元非线性回归模型、相关性分析、亮度差异、暗原色先验、时间序列预测、灰色预测、ARIMA时间序列模型

1．问题重述

1.1 问题背景

能见度是气象、公路行车、飞机飞行中常见指标，单位通常是米。一般意义下的能见度，是指视力正常的人能将目标物从背景中识别出来的最大距离。在科学研究中，能见度是大气环境研究领域的重要课题，是反应大气透明度的一个指标[1]。气象上所讲的大气能见度指白天情境下，正常视力所能看到并且辨认出的以地平线附近天空为背景的大小合适的黑色目标物所对应的最远的距离[2][5]。

能见度对高速公路行车安全非常重要，雾天能见度低对交通运输会产生很大影响，尤其对高速公路运输和航空服务等，在这样的情况下，人们的生命财产安全面临巨大的威胁。据相关部门调查研究发现，发生的交通事故中，由于雾天能见度降低造成严重后果的交通事故数量远远超过其他恶劣天气，且死伤人数占了交通事故总人数的1/3[6]。虽然在公路运输和航空服务中通行标准有所不同，但不可否认，气象条件的影响会给相应的管理和协调带来极大的困扰。因此，能见度预测一直是高速公路部门和航空公司十分重视的问题。

自从 20 世纪初期 Koschmieder 提出 Koschmieder 定律，能见度研究进展逐渐深入。目前，我国高速路网已逐步形成，若大量使用激光能见度仪对全国高速路网进行全覆盖将耗资巨大，同时激光能见度仪还存在对团雾检测精度不高，探测的范围很小，维护成本高等不足。近年来，基于视频的路况（跑道）能见度检测方法受到人们的关注，它在某种程度上克服了激光能见度仪的不足。视频能见度检测方法是将大气光学分析与图像处理及人工智能技术结合，通过对视频图像分析处理，建立视频图像与真实场景之间的关系，再根据图像特征的变化，间接计算出能见度数值。这样做有以下几个优点：第一，通过使用道路监控数据提高了监控相机或者是影像采集设备的使用价值；第二，通过监控视频可以方便地、自动地、长期地、连续地、低成本地获取外界环境数据；第三，这种连续无间断的时间序列数据可以对气象变化进行很好地拟合和预测。然而，现有的基于视频图像的能见度检测方法，由于是间接计算，很难准确地估算能见度。特别地，这些方法中大多数只是选取少量视频、截取图像中的某些固有特征基于 Koschmieder 定律进行估计，并没有充分利用视频的连续信息，估计精度不高，存在较大的改进空间。因此，基于视频数据的能见度检测的研究是一个热点也是难点。

事实上，大雾的形成和消散有其自身的规律，通常与近地层的气象因素有关。而视频资料包含了丰富的信息，特别是涵盖了大雾的变化过程信息。充分利用这些信息，不仅可以提高能见度估计精度，也可以对大雾的消散进行预测，这将对国家、相关部门和人民的生命安全具有重要意义。

1.2 问题提出

本文以上述分析为背景，主要考虑一下几个问题：

（1）根据能见度与地面气象观测站（包括温度、湿度和风速等）的数据，揭示气象因子与能见度之间的数学关系。对各气象因子和能见度的变化规律进行分析，包括变化率和时序性；然后评估各气象因子与能见度之间的相关关系，进而对各因子与能见度进行统计描述，分析之间的函数关系，以此构建出不通气候特征下的多元回归方程；

（2）根据已知的机场视频数据和能见度数据，建立出基于视频数据的能见度估计深度学习模型并进行相应的精度评价。通过对图像特征的提取分析各特征与能见度之间的相关关系，确定深度学习的特征参数向量，并以此构建深度学习模型。

（3）根据提供的视频监控数据，建立通过视频数据估计能见度的算法；由于激光能见度仪不能大规模使用，而高速公路的监控视频数据作为一项简单易得的数据，因此建立不依赖于能见度仪观测数据的能见度观测方法具有十分重要的现实意义；

        （4） 基于前述模型和算法分析能见度随时间变化规律，构建大雾变化的预测模型并根据现实情况下的一些因素对未来一定时间段内的变化曲线做出一定的说明。

2．模型假设

2.1 模型的基本

          在本文的研究中心，采用基于暗通道先验理论所需要进行的问题假设：
          （1）假设研究的道路是光滑的水平道路，即不存在减速装置和道路起伏；
          （2）假设相机拍摄不同时间视频或图像时电子元件均正常；
          （3）假设实际生活中大气本身并不是十分干净，即便是晴天图像中观察远距离物体仍然会存在少量雾霾；

2.2 模型的符号说明

表2-1 模型符号说明

序号	符号	说明
1	s,meas	图像均方差
2	g grad	图像平均梯度
3	c cons	图像对比度
4	µ	图像均值
5	reg	图像子窗格
6	VIS	能见度
7	VIS  ˆ	能见度线性拟合值
8	e	拟合残差
9	ˆe	修正拟合残差
10	vs	风速
11	pa	压强
12	t	温度
13	u	相对湿度

3．技术路线

基于多因素的非线性能见度估计模型

能见度时序变化特

影响能见度变化的气象因子分析

各因子时间序列

征分析

平均 跑道 视程 描述 RVR

平均 气象 光学 视程 描述 MOR

温度

露点温度

相对湿度

本站气压

飞机着陆点气压

修正海平面气压

平均垂直风速

灯光数据

平均风速

平均风向

位置

能见度

单因子选择

统计描述分析

单因子影响程度 分析

能见度与气象因子非线性建模

基于遗传算法的最优参数求解

能见度与气象因子多变量非线性

图3-1 问题一技术流程图

图3-2 问题二技术流程图

4．问题一建模与求解

4.1 问题分析与思路

在问题一中，要求根据题目提供的数据描述能见度与地面气象观测（温度、湿度和风速）的关系，并给出具体的关系式。影响能见度的重要因素是雾和霾，其形成与气象因子具有重要的相关关系，为从定量何定性两个角度去评估和分析二者之间的内在关联，本章尝试利用统计学方法分析某机场气候因子与能见度的特征变化，在此基础上结合数学回归模型构建二者之间的关系表达式，具体做法如下：

首先，对题目提供的AMOS观测数据进行预处理，按照合适的时间间隔关联能见度观测数据与气候数据，然后结合文献调查和相关专业知识对数据中的各项指标进行统计归纳整理，并利用描述性统计方法分析各项因子的特征和变化规律，以去除冗余变量和挖掘各因子对能见度的影响信息分析其间的内在联系，常用的方法有分类统计、倾向变化率、相关系数和可视化分析；然后基于以上特征信息利用非线性回归方程模拟单个气象因子与能见度之间的统计学关系，最后以该统计学关系为基础构建所有气象因子与能见度之间的多元回归模型，利用遗传算法进行最优参数求解从而确定能见度和各气象因子之间的非线性回归方程。

基于多因素的非线性能见度估计模型

能见度时序变化特

影响能见度变化的气象因子分析

各因子时间序列

征分析

平均 跑道 视程 描述 RVR

平均 气象 光学 视程 描述 MOR

温度

露点温度

相对湿度

本站气压

飞机着陆点气压

修正海平面气压

平均垂直风速

灯光数据

平均风速

平均风向

位置

能见度

单因子选择

统计描述分析

单因子影响程度 分析

能见度与气象因子非线性建模

基于遗传算法的最优参数求解

能见度与气象因子多变量非线性

模型

 

图4-1 问题一技术流程图

4.2 数据预处理

4.2.1 信息剔除与数据化

由于气象观测信息很多，形式也多种多样，因而首先需要筛选题目数据中与雾的形成和消散有关的信息，将不必要的数据剔除。共提取出与能见度相关的气象因子如表 4-1 所示。

表4-1 提取因子说明

变量	说明
RVR_1A MOR_1A PAINS(HPA) QFE 06 QNH TEMP RH DEWPOINT WS2A WD2A CW2A	1分钟平均跑道视程 1分钟平均气象光学视程，即能见度 本站气压 飞机着陆地区最高气压 修正海平面气压 温度 相对湿度 露点温度 2分钟平均风速 2分钟平均风向 2分钟平均垂直风速

4.2.2 数据标准化

数据中PTU_R06_15.his的时间间隔为1分钟，而VIS_R06_15.his和WIND_R06_15.his的时间间隔为 15 秒，因此需要统一时间尺度将气象因子与能见度数据进行合并。考虑到数据记录的时间为大约为第一天的8点到第二天8点，时间跨度较大，因而为较好的反应一天的能见度和各气象因子变化以及更方便有效地处理分析数据，最终选取以 15 分钟为间隔整合数据作为后续的实验样本。

表4-2 处理后时间间隔为15分钟的部分数据

LOCALDATE	RVR_	MOR	WS2	WD2A	CW2A	PAINS	QNH	QFE	TEMP	RH	DEWP
(BEIJING)	1A	_1A	A			(HPA)		R06			OINT
2019… 8:00:00	3000	7000	3.78	86	1.78	1023.3	1024.8	1023.2	8.6	90	7.06
2019… 8:15:00	3000	8000	4.45	94	2.61	1023.3	1024.8	1023.2	9.5	86	7.29
2019… 8:30:00	3000	7000	4.55	89	2.34	1023.2	1024.7	1023.1	9.9	84	7.34
2019… 8:45:00	3000	7000	3.73	93	2.14	1023.3	1024.8	1023.2	10.6	82	7.68
2019… 9:00:00	3000	7000	4.12	108	3.15	1023.2	1024.7	1023.1	11.2	80	7.9
2019… 9:15:00	3000	6000	5.7	113	4.67	1023.3	1024.8	1023.2	12.1	78	8.41

4.3 能见度时序变化特征分析
          将两个数据集AMOS20200313和AMOS20191216的能见度数据以15分钟为间隔进行可视化，同时计算其不同时段的倾向变化率，如图4-2所示。

图4-2 能见度时序变化特征图（左图为AMOS20191216，右图AMOS20200313）

从图 4-2 中可以看出，能见度的曲线变化一天的变化趋势主要体现在凌晨附近开始下降，早上8点开始逐渐上升。平均跑道视程曲线的突变时段主要体现在凌晨（00:00）到早上7点附近，气象光学视程整天都在波动变化的突变时段主要为凌晨和早上8点附近。

4.4 气象因子和能见度的相关性分析
          题目所给数据中包含了多种气象数据，为了更加科学地理解和阐述这些变量之间的关系，对能见度和所有气象因子进行 Pearson、Kendall 和 Spearman 检验，检验结果如下图4-3所示。

2019年Pearson相关系数图

2019年Kendall相关系数图

2019年Spearman相关系数图

 

2020年Pearson相关系数图

2020年Kendall相关系数图

2020年Spearman相关系数图

图4-3 各变量的相关性分析

          从图 4-3 中可以看出，气象光学视程和跑道视程之间具有很强的正相关性，查阅相关文献[3][4]可知，气象光学视程是气象学上的一个广义和普遍的概念，即我们常说的能见度

代码

%%   数据与处理

clear;clc;

[datanum,txt] = xlsread("dataAMOS20191216");

datatxt = txt(3:end,2);

newtxt = {};

newNum = [];

for i = 1:size(datatxt,1)

  time = datatxt(i);

  timechar = cell2mat(time);
  timeMilisecond = timechar(length(timechar)-1:length(timechar));%          %1秒间隔   standard = '00';

  if timeMilisecond == standard    %1秒间隔

  newtxt(size(newtxt,1)+1,1) = datatxt(i);

  newNum(size(newNum,1)+1,:) = datanum(i,:);

  end

end

%% 处理15秒一个间隔

clear;clc;
[datanum,txt] = xlsread("dataAMOS20191216",'处理后1分钟间隔'); datatxt = txt(3:end,2);

newtxt = {};

newNum = [];

for i = 1:size(datatxt,1)

  time = datatxt(i);

  timechar = cell2mat(time);
  timeMilisecond = timechar(length(timechar)-4:length(timechar)-3);   if mod(str2num(timeMilisecond),10) ==0 %15秒或30一个间隔   newtxt(size(newtxt,1)+1,1) = datatxt(i);

  newNum(size(newNum,1)+1,:) = datanum(i,:);

  end

end

%遗传算法

clear;clc;
% 目标函数y = a0 + a1*x1^2 + a2*x1 + a3*x2^3 + a4*x2^2 + a5*x2 + a6*x3^4 + a7*x3^3 + a8*x3^2 + a9*x3 + a10*x4^2 + a11*x4
%两边求对数得到化为线性：lny = lna0 + ... + lna11 + 3lnx1 + 6lnx2 + 10lnx3 + 3lnx4 %在令  new = lny; new1 = lnx1; new2 = lnx2; new3 =lnx3; new4 = lnx4;

%得到方程  new = b0 + 3*new1 + 6*new2 + 10*new3 + 3*new4; global    y;

global    x;
adress = {'dataAMOS20200313.xlsx','dataAMOS20191216.xlsx'}; for j = 1:2
data = xlsread(cell2mat(adress(j)),"处理后15分钟间隔");
y = data(:,4)'; %能见度

x1 = data(:,6); %风速

x2 = data(:,9); %本站气压

x3 = data(:,12); %温度

x4 = data(:,13); %相对湿度

x = [x1 x2 x3 x4]';

%定义待调用的两个子程序seaircal()和seair_constraint() ObjectiveFunction = @seaircal;

nvars = 12;

%循环ga()遗传算法

  for i =1:200
            [d,fval] = ga(ObjectiveFunction,nvars,[],[],[],[],[],[],[]); %x是待估变量，fval是返回

的自适应值

  result(i,:) = [d,-fval];

  end

  %将计算结果导出为*.xls 格式的电子表格文档

result1 = sortrows(result,13,'descend');

xlswrite(cell2mat(adress(j)),result1,'遗传算法估计参数');

end
%% 程序功能：为遗传算法提供自适应函数（SEAIIDR模型），返回simerr作为评价指标，即每次进化会选出最小的simerr值

function simerr = seaircal(a)

global x;

global y;

for j = 1:size(x,2)
         yichuan(j) = a(1) + a(2)*x(1,j)^2 + a(3)*x(1,j) + a(4)*x(2,j)^3 + a(5)*x(2,j)^2 + a(6)*x(2,j) + a(7)*x(3,j)^4 + a(8)*x(3,j)^3 + a(9)*x(3,j)^2 + a(10)*x(3,j) + a(11)*x(4,j)^2 + a(12)*x(4,j); end

%% 计算R2值
%计算误差评定指数：模型参数条件下SARS 传播过程与真实过程的拟合度值 simerr = -FITNESS(y',yichuan'); %因为ga()函数默认输出

问题二

% 调试
% Input_path = 'F:\永无极限小组\E20104590100\问题二\图片\裁剪后\';
% namelist = dir(strcat(Input_path,'*.jpg'));
% name=namelist(4).name;    %namelist(i).name; %这里获得的只是该路径下的文件名

% I=imread(strcat(Input_path, name)); %图片完整的路径名

% [meansquare,meanG,tidu,contrast, regionI, regionG, const] = featureExtrude(I,8,8);

%  function  [meansquare,meanG,tidu,contrast,  regionI,  regionG,  const]  =  featureExtrude(image,

ch, cw)

function feature = featureExtrude(image, ch, cw) %%ch和cw是子窗体的尺寸

%获取图像的rgb通道矩阵，并将其转化为double可计算类型      R=image(:,:,1);

  G=image(:,:,2);

  B=image(:,:,3);

  R0=double(R);

  G0=double(G);

  B0=double(B);

  [rs,cs] = size(R0);% 行数列数

  %将rgb颜色空间转换到HSI颜色空间.I为亮度特征值，为了计算均方差   I = 1/3*(R0 + G0 + B0);

  %将rgb颜色空间转换到灰度图像，为了计算梯度   G = double(rgb2gray(image)); %按照0.3 0.69 0.11（rgb）

%*************分块，即按照一定窗口大小计算各个特征***************************   % ch为行间隔  cw为列间隔 ，控制分块图像的大小，如现在的8*8
  % numr为间隔块个数  numc为间隔块个数

  numr = round(rs/ch); %计算总共多少行

  numc = round(cs/cw); %计算总共多少列

  allregion = numr*numc; %总共多少区域

  % 区域块分割

  t1 = (0:numr-1)*ch + 1; t2 = (1:numr)*ch;

  t3 = (0:numc-1)*cw + 1; t4 = (1:numc)*cw;
  count = 1; %用来记录一幅影像生成了多少个窗口，最终结果应该等于allregion+1   %获得一幅影像中的各窗口矩阵，即分块矩阵

  for i=1:numr

  for k=1:numc

  %从t1行至t2的子块

  x=t1(i):t2(i);

  y=t3(k):t4(k);
  regionI(:,:,count) = I(x(1):x(length(x)),y(1):y(length(y)));   regionG(:,:,count) = G(x(1):x(length(x)),y(1):y(length(y)));   count = count + 1;

%                          rectangle ('Position', [t3(j) t1(i) length(x) length(y)], ...,

%                                  'EdgeColor', 'r', 'LineWidth', 1); %绘制矩形分割格子，红色      end

         end
%****************分块，即按照一定窗口大小计算各个特征************************

%计算亮度均方差(是一个窗格的数值)

  for j = 1: allregion

  f = regionI(:,:,j);

  u = mean(mean(f));

  s = 0; %每一个格子的方差

  for i=1:cw

  for k=1:ch

  s = s+(f(i,k)-u)^2;

  end

            end
            meansquare(j) = sqrt(s/(cw*ch)); %一幅影像一个子窗口的一个均方差特征   end   

%计算灰度均值

  for j = 1: allregion

            f = regionG(:,:,j);
            meanG(j) = mean(mean(f));%一幅影像一个子窗口的一个均值特征             end

%计算平均梯度（每个像素的值是梯度，对于窗口还需要求窗口内的所有像素均值,即我们

所做的这个tidu_j/((cw-1)*(ch-1))）

  for j = 1: allregion

  f = regionG(:,:,j);

  tidu_j = 0; %每一个格子的平均梯度

  for i=1:cw-1

  for k=1:ch-1
            tidu_j=(tidu_j)+sqrt( ((f(i,k)-f(i+1,k))^2+(f(i,k)-f(i,k+1))^2)/2 );   end

            end
            tidu(j) = tidu_j/((cw-1)*(ch-1)); %一幅影像一个子窗口的一个均方差特征   end

%计算灰度对比度（每个像素的值是对比度，对于窗口还需要求窗口内的所有像素均值）       for j = 1: allregion

  f = regionG(:,:,j);

  gray_j = []; %每一个格子的对比度

  d = 1;

  for i=2:cw-1

  for k=2:ch-1

  gray_j(1) = abs(f(i,k)-f(i-1,k-1))/max(f(i,k), f(i-1,k-1));

  gray_j(2) = abs(f(i,k)-f(i-1,k))/max(f(i,k), f(i-1,k));

  gray_j(3) = abs(f(i,k)-f(i-1,k+1))/max(f(i,k), f(i-1,k+1));

  gray_j(4) = abs(f(i,k)-f(i,k-1))/max(f(i,k), f(i,k-1));

  gray_j(5) = abs(f(i,k)-f(i,k+1))/max(f(i,k), f(i,k+1));

  gray_j(6) = abs(f(i,k)-f(i+1,k-1))/max(f(i,k), f(i+1,k-1));