2023认证杯A题：太阳黑子预测|数学中国数学建模国际赛（小美赛） |数学建模完整代码+建模过程全解全析_arima模型分析太阳黑子

当大家面临着复杂的数学建模问题时，你是否曾经感到茫然无措？作为2022年美国大学生数学建模比赛的O奖得主，我为大家提供了一套优秀的解题思路，让你轻松应对各种难题。
让我们来看看认证杯的A题！

完整内容可以在文章末尾领取！

Task 1: 请预测当前太阳周期和下一个太阳周期的开始和结束时间；

时间序列分析：

1. 趋势（Trend）： 表示数据长期上升或下降的方向。趋势可以通过拟合多项式或其他非线性函数来建模。在太阳活动中，趋势可能表示太阳黑子数量随时间的总体增加或减少趋势。

2. 季节性（Seasonality）： 表示数据在特定时间范围内的周期性变化。太阳活动通常以约11年的周期性变化为主。

3. 残差（Residuals）： 指的是模型不能解释的数据部分，可能包含随机噪声或其他未考虑的因素。

ARIMA 模型：

ARIMA（差分整合移动平均自回归）是一种常用于时间序列预测的模型。ARIMA模型包括三个主要组件：

1. 自回归（AR）： 表示当前值与过去值的关系，即当前值的预测取决于前一时刻的值。

2. 差分整合（I）： 表示对时间序列进行差分以消除趋势或季节性。

3. 移动平均（MA）： 表示当前值与过去的误差项的关系，即当前值的预测取决于过去的误差。

ARIMA 模型的一般形式为 ARIMA(p, d, q)，其中 p、d、q 分别为自回归、差分整合和移动平均的阶数。

公式示例：

ARIMA 模型的基本预测公式可以表示为：

$$Y_t=c+{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+\dots +{\varphi }_p{Y}_{t-p}+{ϵ}_t-{\theta }_1{ϵ}_{t-1}-{\theta }_2{ϵ}_{t-2}-\dots -{\theta }_q{ϵ}_{t-q}$$

其中：

• $Y_t$ 是时间 t 的观测值。
• $c$ 是常数项。
• ${\varphi }_1,{\varphi }_2,\dots ,{\varphi }_p$ 是自回归系数。
• ${ϵ}_t,{ϵ}_{t-1},\dots ,{ϵ}_{t-q}$ 是残差。

1. 自回归（AR）部分：

   $$Y_t=c+{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+\dots +{\varphi }_p{Y}_{t-p}+{ϵ}_t$$

   这一部分表示当前值 $Y_t$ 与过去 $p$ 个时间步的值 $Y_{t-1},Y_{t-2},\dots ,Y_{t-p}$ 之间的关系。系数 ${\varphi }_1,{\varphi }_2,\dots ,{\varphi }_p$ 表示对应时间步的权重。

2. 差分整合（I）部分：

   差分整合部分表示对时间序列进行 $d$ 阶差分，以消除趋势或季节性。差分整合的目标是使时间序列平稳。

   $$Y_t^{\prime }=Y_t-Y_{t-1}$$

   重复差分 $d$ 次，直到达到平稳性。

3. 移动平均（MA）部分：

   $$Y_t=c+{ϵ}_t-{\theta }_1{ϵ}_{t-1}-{\theta }_2{ϵ}_{t-2}-\dots -{\theta }_q{ϵ}_{t-q}$$

   这一部分表示当前值 $Y_t$ 与过去 $q$ 个时间步的误差项 ${ϵ}_{t-1},{ϵ}_{t-2},\dots ,{ϵ}_{t-q}$ 之间的关系。

使用时间序列分析中的趋势分析方法。其中利用移动平均和趋势分析来估计太阳黑子的趋势。

移动平均法

移动平均法是一种平滑时间序列的方法，可以消除一些噪声，突出趋势。有两种常见的移动平均：简单移动平均（SMA）和指数加权移动平均（EWMA）。

1. 简单移动平均（SMA）：

   $$M{A}_t=\frac{1}{n}\sum _{i=1}^n{Y}_{t-i}$$

   其中 $M{A}_t$ 是时间 $t$ 的移动平均值，$n$ 是选定的时间窗口大小。

2. 指数加权移动平均（EWMA）：

   $$EM{A}_t=\alpha \cdot Y_t+(1-\alpha )\cdot EM{A}_{t-1}$$

   其中 $EM{A}_t$ 是时间 $t$ 的指数加权移动平均值，$\alpha$ 是平滑因子$（0<\alpha <1）$。

趋势分析

趋势分析旨在捕捉数据中的长期趋势。可以通过拟合多项式或其他函数来实现。假设我们采用线性趋势：

$$Y_t={\beta }_0+{\beta }_1\cdot t+{ϵ}_t$$

其中 $Y_t$ 是时间 $t$ 的观测值，${\beta }_0$ 是截距，${\beta }_1$ 是斜率，${ϵ}_t$ 是误差项。

整合趋势和移动平均

可以通过以下步骤来整合趋势和移动平均来预测太阳黑子的趋势：

1. 计算移动平均： 使用简单移动平均或指数加权移动平均平滑太阳黑子的时间序列数据。

2. 趋势分析： 利用线性趋势分析，拟合趋势方程。

3. 整合： 将移动平均和趋势相加，得到整合后的趋势估计。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# 生成模拟太阳黑子数据
np.random.seed(42)
t = np.arange(1, 101)
sunspots = 10 * np.sin(0.2 * t) + np.random.normal(scale=2, size=len(t))

# 创建DataFrame
data = pd.DataFrame({'Time': t, 'Sunspots': sunspots})

# 计算简单移动平均
window_size = 5
data['SMA'] = data['Sunspots'].rolling(window=window_size).mean()

# 趋势分析
X = data['Time'].values.reshape(-1, 1)
y = data['Sunspots'].values

model = LinearRegression()
model.fit(X, y)
data['Trend'] = model.predict(X)

# 整合趋势和移动平均
data['Integrated_Trend'] = data['SMA'] + data['Trend']

# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(data['Time'], data['Sunspots'], label='Original Data')
plt.plot(data['Time'], data['SMA'], label=f'SMA ({window_size} periods)')
plt.plot(data['Time'], data['Trend'], label='Trend')
plt.plot(data['Time'], data['Integrated_Trend'], label='Integrated Trend')

plt.xlabel('Time')
plt.ylabel('Sunspots')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

Task 2: 请预测当前太阳周期和下一个太阳周期的开始和结束时间；

理论基础：

1. 太阳活动周期： 太阳活动表现为周期性的变化，主要是由于太阳磁场的活动导致的。太阳活动的周期性变化约为11年，被称为太阳活动周期或太阳周期。这个周期性变化主要是由太阳黑子的数量和太阳黑子区域的活动引起的。

2. 太阳黑子与太阳活动关系： 太阳黑子是太阳表面的一种现象，它们标志着强烈的磁活动。太阳黑子的数量和活动水平通常与整个太阳活动周期相关联，即在太阳周期内，太阳黑子的数量经历峰值和谷值的变化。

3. 峰值时期： 太阳周期中的峰值时期被称为太阳最大值。在太阳最大值时，太阳黑子的数量和太阳黑子区域的活动都达到峰值。这通常是太阳活动的高峰期，伴随着较强的太阳射击和其他太阳活动。

4. 太阳最大值的预测： 预测太阳最大值的发生时间和持续时间是一个复杂的问题，因为它涉及到理解太阳活动的驱动机制以及对太阳周期的整体变化的建模。传统上，科学家们使用太阳黑子的观测数据，结合数学模型，例如傅里叶分析、差分方程、时间序列分析等，来尝试预测太阳最大值。

5.峰值的持续时间： 太阳最大值的持续时间通常取决于太阳黑子的活动水平和持续时间。在太阳最大值期间，太阳黑子活动较为频繁和强烈。持续时间的估计可能涉及到对太阳黑子数据的统计分析和模型拟合。

预测太阳活动的确切时刻和强度是一个非常复杂的问题，因为太阳活动受多种影响，包括太阳磁场的变化和其他不确定性因素。目前，没有单一的理论公式可以准确描述太阳活动周期的最大值时刻和持续时间。然而，一些传统的方法和模型被用来尝试进行这样的预测。

经典方法：

1. Wolf数（Wolf Number）：

   Wolf数是一种用于衡量太阳黑子活动水平的指标，可以用来估计太阳活动周期的强度。它通常由以下公式定义：

   $$W=10\cdot G+N$$

   其中 $ G $ 是给定时期内的太阳黑子群的数量，$ N $ 是给定时期内的太阳黑子的总数。

2. 斯皮格勒公式（Spörer’s Law）：

   斯皮格勒公式描述了太阳活动周期的变化，其基本形式为：

   $$D=A+B\cdot C^M$$

   其中 $ D $ 是太阳活动周期的持续时间，$ A $、$ B $、$ C $ 和 $ M $ 是与具体时期相关的参数。

正弦函数是一种简单而常用的周期性函数，可以用来描述太阳活动的周期性变化。在这里，我们将使用正弦函数来进行太阳活动周期的分析：

$$Y(t)=A\sin (2\pi ft+\varphi )+C$$

其中：

• $ Y(t) $ 是太阳黑子的数量或其他太阳活动指标。
• $ A $ 是振幅，表示活动的强度。
• $ f $ 是频率，与太阳周期的倒数相关。
• $ t $ 是时间。
• $ \phi $ 是相位，表示正弦波的偏移。
• $ C $ 是偏移项，表示整体的水平偏移。

Python 代码示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 模拟太阳活动周期数据
np.random.seed(42)
t = np.arange(1, 301)
sunspots = 10 * np.sin(0.1 * t) + np.random.normal(scale=2, size=len(t))

# 定义正弦函数模型
def sine_function(t, A, f, phi, C):
    return A * np.sin(2 * np.pi * f * t + phi) + C

# 使用 curve_fit 拟合正弦函数模型
params, covariance = curve_fit(sine_function, t, sunspots)

# 提取拟合后的参数
A_fit, f_fit, phi_fit, C_fit = params

# 生成拟合后的曲线
fit_curve = sine_function(t, A_fit, f_fit, phi_fit, C_fit)

# 绘制原始数据和拟合后的曲线
plt.figure(figsize=(10, 6))
plt.plot(t, sunspots, label='Original Data')
plt.plot(t, fit_curve, label='Fit Curve')

plt.xlabel('Time')
plt.ylabel('Sunspots')
plt.legend()
plt.show()

# 输出拟合后的参数
print(f"振幅 (A): {A_fit}")
print(f"频率 (f): {f_fit}")
print(f"相位 (phi): {phi_fit}")
print(f"偏移项 (C): {C_fit}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

针对太阳周期的最大值，可以使用峰值检测方法，找到太阳活动周期中的峰值时刻。并观察活动水平的持续时间。

import numpy as np
import pandas as pd
from scipy.signal import find_peaks

# 生成模拟太阳黑子数据
np.random.seed(42)
t = np.arange(1, 301)
sunspots = 10 * np.sin(0.1 * t) + np.random.normal(scale=2, size=len(t))

# 创建DataFrame
data = pd.DataFrame({'Time': t, 'Sunspots': sunspots})

# 寻找峰值
peaks, _ = find_peaks(data['Sunspots'], height=0)

# 提取峰值对应的时间点
peak_times = data['Time'].iloc[peaks]

# 计算持续时间
if len(peak_times) >= 2:
    max_intensity_start = peak_times.iloc[0]
    max_intensity_end = peak_times.iloc[-1]
    max_intensity_duration = max_intensity_end - max_intensity_start
else:
    max_intensity_start = max_intensity_end = max_intensity_duration = np.nan

print(f"预测太阳周期太阳最大值的发生时间：从 {max_intensity_start} 到 {max_intensity_end}")
print(f"预测太阳周期太阳最大值的持续时间：{max_intensity_duration} 单位时间")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

Task 3: 预测当前太阳周期和下一个太阳周期的太阳黑子数量和面积，并在您的论文中解释模型的可靠性。

太阳周期信息： 使用历史数据分析，确定太阳周期的长度和变化。

历史太阳黑子数量： 考虑前几个周期的太阳黑子数量作为特征。

我们使用长短时记忆网络（LSTM）进行时间序列预测时，模型的基本公式可以描述为以下过程：

1. 输入层：

   输入层接受时间序列数据，其形状为 $(batch\_size,time\_steps,features)$。

2. LSTM 层：

   LSTM 层是模型的核心部分，其内部包含许多门（gates），用于控制信息的流动。LSTM 的更新过程如下：

   • 输入门（Input Gate）：
     $$i_t=\sigma (W_{ii}{x}_t+b_{ii}+W_{hi}{h}_{t-1}+b_{hi})$$

   • 遗忘门（Forget Gate）：
     $$f_t=\sigma (W_{if}{x}_t+b_{if}+W_{hf}{h}_{t-1}+b_{hf})$$

   • 细胞状态（Cell State）更新：
     $$g_t=\tanh (W_{ig}{x}_t+b_{ig}+W_{hg}{h}_{t-1}+b_{hg})$$
     $$c_t=f_t\cdot c_{t-1}+i_t\cdot g_t$$

   • 输出门（Output Gate）：
     $$o_t=\sigma (W_{io}{x}_t+b_{io}+W_{ho}{h}_{t-1}+b_{ho})$$

   • 隐藏状态更新：
     $$h_t=o_t\cdot \tanh (c_t)$$

   其中，$\sigma$ 表示 sigmoid 激活函数，$\tanh$ 表示双曲正切激活函数，$W$ 和 $b$ 分别是权重和偏置。

3. 输出层：

   输出层产生模型的最终预测结果。在时间序列预测中，可以是单个值或多个值，具体取决于任务的性质。

4. 损失函数：

   损失函数用于度量模型预测值与实际值之间的差异。在回归任务中，常用均方误差（Mean Squared Error，MSE）作为损失函数。

整个 LSTM 模型的训练过程涉及通过反向传播算法调整权重和偏置，使损失函数最小化。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
from keras.models import Sequential
from keras.layers import LSTM, Dense

# 生成模拟太阳黑子数据
np.random.seed(42)
t = np.arange(1, 501)
sunspots = 10 * np.sin(0.1 * t) + np.random.normal(scale=2, size=len(t))

# 归一化数据
scaler = MinMaxScaler(feature_range=(0, 1))
sunspots_normalized = scaler.fit_transform(sunspots.reshape(-1, 1))

# 创建时间窗口数据集
def create_dataset(dataset, time_steps=1):
    dataX, dataY = [], []
    for i in range(len(dataset) - time_steps):
        a = dataset[i:(i + time_steps), 0]
        dataX.append(a)
        dataY.append(dataset[i + time_steps, 0])
    return np.array(dataX), np.array(dataY)

# 设置时间窗口大小
time_steps = 10

# 创建时间窗口数据集
X, y = create_dataset(sunspots_normalized, time_steps)

# 将数据整形为 LSTM 期望的格式 [样本数, 时间步数, 特征数]
X = np.reshape(X, (X.shape[0], X.shape[1], 1))

# 构建 LSTM 模型
model = Sequential()
model.add(LSTM(units=50, return_sequences=True, input_shape=(X.shape[1], 1)))
model.add(LSTM(units=50))
model.add(Dense(units=1))
model.compile(optimizer='adam', loss='mean_squared_error')

# 训练模型
model.fit(X, y, epochs=100, batch_size=1, verbose=2)

# 使用训练好的模型进行预测
train_predict = model.predict(X)
train_predict = scaler.inverse_transform(train_predict)
y_actual = scaler.inverse_transform(y.reshape(-1, 1))

# 计算均方根误差
rmse = np.sqrt(mean_squared_error(y_actual, train_predict))
print(f"均方根误差 (RMSE): {rmse}")

# 可视化预测结果
plt.figure(figsize=(15, 6))
plt.plot(t[time_steps:], sunspots[time_steps:], label='实际值')
plt.plot(t[time_steps:], train_predict, label='预测值')
plt.xlabel('时间步数')
plt.ylabel('太阳黑子数量')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

可行性分析

在数据预处理阶段，我们对太阳黑子数量和面积的历史数据进行了可行性分析，并采取了一系列合理的处理方法。可行性结果的重新组织：

数据可行性结果：

1. 缺失值和异常值处理：

   • 通过采用线性插值和异常值修正的策略，我们成功处理了历史太阳黑子数据中的缺失值和异常值。这确保了数据的完整性和可用性。

2. 时间序列平滑：

   • 应用了移动平均方法对太阳黑子数量的时间序列进行了平滑处理。这有助于降低数据中的噪声，使其更适合作为模型的输入。

数据处理公式：

# 缺失值处理 - 线性插值
def interpolate_missing_values(data):
    interpolated_data = data.interpolate(method='linear')
    return interpolated_data

# 异常值处理 - 修正为阈值
def handle_outliers(data, threshold):
    data[data > threshold] = threshold
    return data

# 时间序列平滑 - 移动平均
def smooth_time_series(data, window_size):
    smoothed_data = data.rolling(window=window_size).mean()
    return smoothed_data
1
2
3
4
5
6
7
8
9
10
11
12
13
14

评估：

通过对比处理前后数据的统计特征、可视化结果以及在一小部分数据上的模型预测效果，我们对处理后数据的质量和对模型的影响进行了评估。这一系列处理方法在维护数据质量的同时，为后续的预测模型提供了更为可靠的输入。

这样的数据预处理工作为进一步建立和训练预测模型奠定了基础，确保模型在真实世界的应用中能够取得更好的性能。

更多内容具体可以看看我的下方名片！里面包含有认证杯一手资料与分析！
另外在赛中，我们也会陪大家一起解析认证杯的一些方向
关注 CS数模 团队，数模不迷路～