自回归模型是统计学上常用的一类时间序列模型，用一个变量𝑦𝑡的历史信息来预测自己
有外部数据的非线性自回归模型是自回归模型的扩展，在每个时刻 t 都有一个外部输入𝑥𝑡,产生一个输出𝑦𝑡,有外部输入的非线性回归模型通过一个延时器记录最近𝐾𝑥次的外部输入和最近𝐾𝑦次的输出，第 t 个时刻的输出𝑦𝑡为
其中ƒ(·)表示非线性函数，可以使一个前馈网络, 𝐾𝑥和𝐾𝑦为超参数
简单循环网络只有一个隐藏层的神经网络。在一个两层的前馈神经网络中，连接存在相邻的
层与层之间，隐藏层的节点之间是无连接的。而简单循环网络增加了从隐藏层到隐藏层的反
馈连接。
令向量𝑥𝑡 ∈ 𝑅𝑀表示在时刻 t时网络的输入，ℎ𝑡 ∈ 𝑅𝐷表示隐藏层状态(即隐藏层神经元活性值)，则ℎ𝑡不仅和当前时刻的输入𝑥𝑡相关,也和上一个时刻的隐藏层状态ℎ𝑡−1相关。简单循环网络在时刻 t 的更新公式为
如果我们把每个时刻的状态都看作前馈神经网络的一层，循环神经网络可以看作在时间维度上权值共享的神经网络。
我们先定义一个完全连接的循环神经网络，其输入为𝑥𝑡，输入为𝑦𝑡
其中 h 为隐状态，ƒ(·)为非线性激活函数，U、W、b 和 V 为网络参数
以2020年华为杯研究生数学建模竞赛中的E题为例：
原题再现：
为了估计不同大雾情况下对应的能见度以及预测大雾的消散，请回答以下问题：
一． 众所周知，雾与近地面的气象因素有关。建立模型描述能见度与地面气象观测（温度、湿度和风速等）之间的关系，并针对题目所提供的数据（机场AMOS观测.zip）导出具体的关系式；
二． 根据题目提供的某机场视频数据（机场视频.zip）和能见度数据（机场AMOS观测.zip），建立基于视频数据的能见度估计深度学习模型，并对估计的能见度进行精度评估；
任务二提供了某机场的视频数据和 AMOS 观测的能见度数据，要求建立基于视频数据的能见度估计深度学习模型，并对其精度进行评估。视频数据是是由一组连续的图像构成的，其每幅图像具有时间属性，也表现出一定的空间范围的信息，属于常见的时空数据。通过视频数据可知，每秒有 25 帧图像，时间跨度是 2020 年 3 月 13 日 0 时 0 分 26 秒开始，当日 11 时 47 分 48 秒结束。视频的观测点是固定的，但是视频的观测视角和视野范围约在 9 时 13 分后存在变化。另外，由于雾气较重、能见度较低，视频部分数据无法体现景物信息，目标观测物体已完全被景物遮挡。需要指出的是，由于观测设备质量、大气颗粒物等影响，视频的部分图像存在色彩、亮度突变的异常情况，同时也存在大量的噪声点。视频数据属于二维数据，而机场 AMOS 数据属于一维数据，时间分辨率为 15 秒，时间跨度是 2020 年 3 月 13 日 0 时 0 分 26 秒开始，次日 7 时 59 分 45 秒结束，两者的时间跨度和分辨率和不同。基于以上的情况分析，需要根据 AMOS 能见度观测时间提取视频中的关键帧，以此将某时刻的图像和能见度进行匹配。
为了建立视频数据和能见度数据之间的深度学习模型，神经网络（Artificial Neural Networks，ANN）是必不可少的工具，其基于梯度下降进行参数寻优，具有较强的学习能力。通常对于二维图像数据，卷积神经网络（Convolutional Neural Networks，CNN）对于图像结构化特征提取具有很大优势，也广泛应用于目标检测、语义分割等场景。但由于本题的图像数据存在上述异常情况，其局部特征对于能见度的反演能力有限，并且图像数据量有限，直接基于图像进行深度学习并不合适。因此，需要对图像数据进行预处理，以提高图像的质量，增强图像的特征，减少噪声等随机信息对模型收敛的影响。基于特征工程的思想，本文先提取每一幅图像的特征信息，再采用全连接神经网络进行学习，这样的深度学习模型具有收敛速度快、计算成本低、可解释性高的优势。
首先对给出的机场视频数据（Fog20200313000026.mp4）进行分析，视频总时长 11 时33 分 28 秒，视频分辨率为 1280×720，图像为三通道图像，帧率为 25 帧/秒，拍摄的时间跨度是 2020 年 3 月 13 日 0 时 0 分 26 秒开始，当日 11 时 47 分 48 秒结束。通过 Adobe Premiere Pro 软件观察数据发现，该视频帧率和左上角水印时间匹配存在些许差异，每秒约为 24~26 帧。为了将视频数据与能见度数据进行时间上的匹配，首先获取能见度和视频的共同观测部分，即从 2020 年 3 月 13 日凌晨 0 时 0 分 30 秒至上午 7 时 59 分 45 秒。其中，在 0 时 47 分至 1 时 1 分间视频数据缺失，另在 1 时 39 分 45 秒、4 时 11 分 15 秒 AMOS数据缺失的情况。在此基础上，根据 AMOS 能见度的时间分辨率，时间（每隔 15 秒）提取出重要的关键帧，最终将连续的视频数据转换为1860幅图像，分辨率仍保留为1280×720。 通过这 1860 幅图像可以发现，在大雾天气下，可能由于视频传感器通道和质量问题，图像出现了紫色、蓝色、黄色等异常情况，如图 1 所示。另外，大雾天气下，水汽和气溶胶颗粒分布在空气中，对夜间灯光会产生不稳定的散射，相邻两幅图像的背景灯光亮度会产生较大变化，如图 2所示
图像预处理（去噪、增强、掩膜）
基于以上分析，首先利用不同滤波器对图像进行去噪，通过比较确定适用于的本题的滤波算法。一般图像的噪声可以分为椒盐噪声和高斯噪声，椒盐噪声就是在图像上随机出现黑色白色的像素；高斯噪声是指它的概率密度函数服从高斯分布（即正态分布）的一类噪声，与椒盐噪声相似，区别在于椒盐噪声是出现在随机位置、噪点深度基本固定的噪声，高斯噪声与其相反，是几乎每个点上出现噪声、噪点深度随机的噪声。本题的数据同时存在两种噪声，且不同图像的噪声分布范围大且不稳定。图像去噪是指去除噪声图像中的噪声，从而还原真实图像。但是，由于噪声、边缘和纹理都是高频分量，在去噪过程中很难区分它们，因此去噪后的图像不可避免地会丢失一些细节。总体而言，从有噪声的图像中恢复有意义的信息以获得高质量的图像是一个重要的问题，分别采用中值滤波、均值滤波和高斯滤波，其去噪结果如图 所示。中值滤波是一种非线性空间滤波器，它将图像滤波器包围的区域中像素的统计排序的中位数值代替中心像素的值；均值滤波器的输出是包含在滤波掩模领域内像素的简单平均值；高斯滤波是一种线性平滑滤波，整幅图像的每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均计算后得到，适用于消除高斯噪声。但是三种方法对于能见度较低的图像并不适用。
非局部均值（Non-Local Means）是近年来提出的一项新的去噪算法，其充分利用了图像中的冗余信息，并且能最大程度地保持图像的细节特征。基本思想是利用整幅图像进行去噪，它以图像块为单位在图像中寻找相似区域，再对这些区域进行加权平均平均，能够较好地滤除图像中的噪声。该算法适用于短时间的图像序列，对本题去噪具有较大的优势。
在与视频成功匹配的结果共 1860 条能见度数据中，可以看到能见度的全部样本数据为 0-50m 区间内占比极高（如图），高达 78.6%，而 1000m 以上能见度的样本只占全部数据的 1.5%。需要特别强调的是，训练集的样本不均衡势必会导致深度学习的预测结果趋于样本多的那类情况，这样的模型虽然拟合精度高，但其实属于一种错误的学习，其并不能有效得将图像与能见度建立联系。由于深度学习的学习能力很强，很容易出现过拟合现象，而样本的不均衡可能掩盖了这种过拟合现象。为了得到可靠的深度学习模型，输入的训练数据至关重要，本文先去除大部分 50m能见度的样本，以保证训练数据的能见度分布较为均衡。但是，将所有 50m 能见度的样本作为测试集来评价模型，是不合理的。因此，在样本均衡的能见度中进行分层抽样，将 70%的样本作为真正的训练集，剩余的 30%样本归入测试集，这样的测试集也包含了不同的能见度样本。
经过图像特征点的提取，训练样本和测试样本的划分后，即获取了深度学习模型的输入数据。传统 BP 神经网络算法及其改进算法都是非完全全连接神经网络算法，具有收敛速度慢，泛化能力差等不足。本文构建的是全连接神经网络模型，为了防止过拟合的情况发生，我们在全连接层之间设置了 Dropout 层，使得在训练过程中使部分连接层停止工作，以避免其对于某些局部特征的过度依赖，这样可以进一步提高模型的泛化能力。经过 200 轮次的训练，得到了全连接神经网络模型在训练集和验证集上的损失函数变化曲线（下图所示），损失函数呈现明显的下降趋势，并随着训练的深入达到一个平衡的状态，可以看到在训练集和验证集上，损失函数的表现具有一致性，说明我们的模型并没有出现明显的过拟合现象。
import cv2
import pandas as pd
import numpy as np
import os
from sklearn import svm
df = pd.read_excel('data.xlsx')
y = df['MOR_1A']
X = []
filenames = sorted((fn for fn in os.listdir('C:\\python处理分析excel\\') if fn.endswith('.jpg')))
filenames.sort(key=lambda x: int(x[0:-4]))
sift = cv2.xfeatures2d.SIFT_create()
for file in filenames[0:1]:
img = cv2.imread('C:\\python处理分析excel\\' + file)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
kp = sift.detect(gray, None)
# 显示关键点
ret = cv2.drawKeypoints(gray, kp, img)
cv2.imshow('ret', ret)
cv2.waitKey(0)
cv2.destroyAllWindows()
for file in filenames:
img = cv2.imread('C:\\python处理分析excel\\+ file)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
# 找出关键点并绘图
kp, des = sift.compute(gray, kp)
des = des.reshape((-1,))
# 使用关键点找出 SIFT 特征向量
X.append(list(des))
X = pd.DataFrame(X)
# 划分训练集和测试集
from keras.preprocessing import sequence
from keras.models import Sequential
from keras.datasets import boston_housing
from keras.layers import Dense, Dropout
from keras.utils import multi_gpu_model
from keras import regularizers # 正则化
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
x_train, x_test, y_train, y_test = train_test_split(X, y,
test_size=0.3, stratify=y,
random_state=2020)
model = Sequential() # 初始化，很重要！
model.add(Dense(units=10, # 输出大小
activation='relu', # 激励函数
input_shape=(x_train.shape[1],) # 输入大小, 也就是列的大小
)
)
model.add(Dropout(0.2)) # 丢弃神经元链接概率
model.add(Dense(units=15,
kernel_regularizer=regularizers.l2(0.01), # 施加在权重上的正则项
activity_regularizer=regularizers.l1(0.01), # 施加在输出上的正则项
activation='relu', # 激励函数
bias_regularizer=keras.regularizers.l1_l2(0.01) # 施加在偏置向量上的正则项
)
)
model.add(Dense(units=1,
activation='linear' # 线性激励函数 回归一般在输出层用这个激励函数
)
)
print(model.summary()) # 打印网络层次结构
model.compile(loss='mse', # 损失均方误差
optimizer='adam', # 优化器
)
history = model.fit(x_train, y_train,
epochs=200, # 迭代次数
batch_size=200, # 每次用来梯度下降的批处理数据大小
verbose=2, # verbose：日志冗长度，int：冗长度，0：不输出训练过程，1：输出训练进度，2：输出每一个 epoch
validation_data=(x_test, y_test) # 验证集
)
fig = plt.figure(figsize=(8, 4))
ax = fig.add_subplot(111)
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.ylabel('损失函数', fontsize=12)
plt.xlabel('迭代次数', fontsize=12)
plt.legend(['训练集', '测试集'], loc='upper right', fontsize=12)
plt.savefig('损失曲线.jpg', dpi=600, bbox_inches='tight')
plt.show()
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error
# 散点图
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(211)
predict_y = model.predict(x_train)
RMSE = np.sqrt(mean_squared_error(y_train, predict_y))
MAE = mean_absolute_error(y_train, predict_y)
R2 = r2_score(y_train, predict_y)
plt.scatter(y_train, predict_y, c='red', s=5)
plt.xlabel('[训练集]气象光学视程(m)观测值', fontsize=14)
plt.ylabel('[训练集]气象光学视程(m)模拟值', fontsize=14)
plt.text(0, 900, 'RMSE={},\nMAE={},\n$R^2$={}'.format(int(RMSE), int(MAE), round(R2, 2fontsize=14)
ax = fig.add_subplot(212)
predict_y = model.predict(x_test)
RMSE = np.sqrt(mean_squared_error(y_test, predict_y))
MAE = mean_absolute_error(y_test, predict_y)
R2 = r2_score(y_test, predict_y)
plt.scatter(y_test, predict_y, c='blue', s=5)
plt.xlabel('[测试集]气象光学视程(m)观测值', fontsize=14)
plt.ylabel('[测试集]气象光学视程(m)模拟值', fontsize=14)
plt.text(0, 900, 'RMSE={},\nMAE={},\n$R^2$={}'.format(int(RMSE), int(MAE), round(R2, 2)), fontsize=14)
plt.savefig('模拟结果.jpg', dpi=600, bbox_inches='tight')
plt.savefig('模拟结果.pdf', bbox_inches='tight')
plt.show()
