2023年西南大学数学建模C题天气预报解题全过程文档及程序

2023年西南大学数学建模

C题 天气预报

原题再现：

  天气现象与人类的生产生活、社会经济、军事活动等方方面面都密切相关，大到国家，小到个人，都受到极端天气的影响。2022年6月，全球陆地地区出现了自1850年代末人类有系统气象记录以来最热的温度，干旱、野火、暴雨和洪水席卷北半球。2022年夏季，我国长江流域发生了1961年有记录以来持续时间最长、强度最大的特大干旱。其中，重庆全市出现持续高温天气，多地持续40℃高温。自2022年8月17日重庆涪陵区首次发生山火灾害，十天时间内南川、江津、大足、铜梁、北碚、巴南、璧山、开州、长寿等地相继发生多起山火。
  天气预报是对未来一定时期内天气（如气温、风、云、降水和天气现象）变化的事先估计和预告。天气预报是气象工作为国民经济和国防建设服务的重要手段。天气预报按预报范围而分，有区域天气预报和局地天气预报。按服务对象的需要而分有日常性的公众天气预报和专业天气预报。按时效的长短而分，有短期、中期和长期天气预报。气象台通过各种渠道及时准确地公开发布天气预报，特别是灾害性天气预报，在保护人民生命财产，促进经济发展等方面发挥着重要作用。近代以来，天气预报一直是以天气学原理为基础的。随着计算技术及探测技术的发展，除常规天气图方法结合数理统计方法制作预报外，又将气象雷达和卫星探测资料应用于预报业务，同时发展了数值预报方法。该方法通过确定大气质量、能量和动量的守恒原理来预报大气的物理过程，显著地提高天气形势预报的质量，从而促进天气预报的客观定量化。
  请通过提供的网址或其它资料收集数据，解决如下问题：
  （1）从机理分析的角度，建立不同天气预报策略（预报范围和时效）的数学模型。
  （2）请收集数据，对2022年夏天重庆市的天气，以及发生过特大山火区域的天气进行预测，分析其发生发展的必然性与偶然性；
  （3）请收集数据，以提交论文的时间为起点，对重庆市沙坪坝区的天气进行预测，给出前3天内、第10-15天、第30-35的预测结果。

问题背景：

  天气现象与人类的生产生活、社会经济、军事活动等方方面面都密切相关，大到国家，小到个人，都受到极端天气的影响。2022年6月，全球陆地地区出现了自1850年代末人类有系统气象记录以来最热的温度，干旱、野火、暴雨和洪水席卷北半球。2022年夏季，我国长江流域发生了1961年有记录以来持续时间最长、强度最大的特大干旱。其中，重庆全市出现持续高温天气，多地持续40℃高温。自2022年8月17日重庆涪陵区首次发生山火灾害，十天时间内南川、江津、大足、铜梁、北碚、巴南、璧山、开州、长寿等地相继发生多起山火。

问题重述：

  请通过提供的网址或其它资料收集数据，解决如下问题：
  （1）根据实际采集数据，从多角度考虑，建立不同天气预报策略（预报范围和时效）的数学模型。
  （2）收集相关当地天气数据，对2022年夏天重庆市的天气，以及发生过特大山火区域的天气进行预测，分析其发生发展的必然性与偶然性；
  （3）请收集数据，以提交论文的时间为起点，对重庆市沙坪坝区的天气进行预测，给出前3天内、第10-15天、第30-35的预测结果。

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

非慈善耶稣

模型的建立与求解整体论文缩略图

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

非慈善耶稣

部分程序代码：

i

mport numpy as np; 
import statsmodels.api as sm
a=np.loadtxt("Pdata12_3.txt")   #加载表中x1,x2,x3,y的11行4列数据
x=a[:,:3]  #提出自变量观测值矩阵
X=sm.add_constant(x)  #增加第一列全部元素为1得到增广矩阵
md=sm.OLS(a[:,3],X).fit()  #构建并拟合模型
b=md.params          #提取所有回归系数
y=md.predict(X)      #求已知自变量值的预测值
print(md.summary())  #输出模型的所有结果
print("相关系数矩阵:\n",np.corrcoef(x.T))
X1=sm.add_constant(a[:,0])
md1=sm.OLS(a[:,2],X1).fit()
print("回归系数为：",md1.params)



灰色预测
from decimal import *


class GM11():
    def __init__(self):
        self.f = None

    def isUsable(self, X0):
        '''判断是否通过光滑检验'''
        X1 = X0.cumsum()
        rho = [X0[i] / X1[i - 1] for i in range(1, len(X0))]
        rho_ratio = [rho[i + 1] / rho[i] for i in range(len(rho) - 1)]
        print("rho:", rho)
        print("rho_ratio:", rho_ratio)
        flag = True
        for i in range(2, len(rho) - 1):
            if rho[i] > 0.5 or rho[i + 1] / rho[i] >= 1:
                flag = False
        if rho[-1] > 0.5:
            flag = False
        if flag:
            print("数据通过光滑校验")
        else:
            print("该数据未通过光滑校验")

        '''判断是否通过级比检验'''
        lambds = [X0[i - 1] / X0[i] for i in range(1, len(X0))]
        X_min = np.e ** (-2 / (len(X0) + 1))
        X_max = np.e ** (2 / (len(X0) + 1))
        for lambd in lambds:
            if lambd < X_min or lambd > X_max:
                print('该数据未通过级比检验')
                return
        print('该数据通过级比检验')

    def train(self, X0):
        X1 = X0.cumsum()
        Z = (np.array([-0.5 * (X1[k - 1] + X1[k]) for k in range(1, len(X1))])).reshape(len(X1) - 1, 1)
        # 数据矩阵A、B
        A = (X0[1:]).reshape(len(Z), 1)
        B = np.hstack((Z, np.ones(len(Z)).reshape(len(Z), 1)))
        # 求灰参数
        a, u = np.linalg.inv(np.matmul(B.T, B)).dot(B.T).dot(A)
        u = Decimal(u[0])
        a = Decimal(a[0])
        print("灰参数a：", a, "，灰参数u：", u)
        self.f = lambda k: (Decimal(X0[0]) - u / a) * np.exp(-a * k) + u / a

    def predict(self, k):
        X1_hat = [float(self.f(k)) for k in range(k)]
        X0_hat = np.diff(X1_hat)
        X0_hat = np.hstack((X1_hat[0], X0_hat))
        return X0_hat

    def evaluate(self, X0_hat, X0):
        '''
        根据后验差比及小误差概率判断预测结果
        :param X0_hat: 预测结果
        :return:
        '''
        S1 = np.std(X0, ddof=1)  # 原始数据样本标准差
        S2 = np.std(X0 - X0_hat, ddof=1)  # 残差数据样本标准差
        C = S2 / S1  # 后验差比
        Pe = np.mean(X0 - X0_hat)
        temp = np.abs((X0 - X0_hat - Pe)) < 0.6745 * S1
        p = np.count_nonzero(temp) / len(X0)  # 计算小误差概率
        print("原数据样本标准差：", S1)
        print("残差样本标准差：", S2)
        print("后验差比：", C)
        print("小误差概率p：", p)


if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import numpy as np

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
    plt.rcParams['axes.unicode_minus'] = False  # 步骤二（解决坐标轴负数的负号显示问题）

    # 原始数据X
    X = np.array(
        [21.2, 22.7, 24.36, 26.22, 28.18, 30.16, 32.34, 34.72, 37.3, 40.34, 44.08, 47.92, 51.96, 56.02, 60.14,
         64.58,
         68.92, 73.36, 78.98, 86.6])
    # 训练集
    X_train = X[:int(len(X) * 0.7)]
    # 测试集
    X_test = X[int(len(X) * 0.7):]

    model = GM11()
    model.isUsable(X_train)  # 判断模型可行性
    model.train(X_train)  # 训练
    Y_pred = model.predict(len(X))  # 预测
    Y_train_pred = Y_pred[:len(X_train)]
    Y_test_pred = Y_pred[len(X_train):]
    score_test = model.evaluate(Y_test_pred, X_test)  # 评估

    # 可视化
    plt.grid()
    plt.plot(np.arange(len(X_train)), X_train, '->')
    plt.plot(np.arange(len(X_train)), Y_train_pred, '-o')
    plt.legend(['降雨实际值', '灰色预测模型预测值'])
    plt.title('训练集')
    plt.show()

    plt.grid()
    plt.plot(np.arange(len(X_test)), X_test, '->')
    plt.plot(np.arange(len(X_test)), Y_test_pred, '-o')
    plt.legend(['负荷实际值', '灰色预测模型预测值'])
    plt.title('测试集')
    plt.show()

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

非慈善耶稣