2023年国赛高教杯数学建模
E题 黄河水沙监测数据分析
原题再现
黄河是中华民族的母亲河。研究黄河水沙通量的变化规律对沿黄流域的环境治理、气候变化和人民生活的影响，以及对优化黄河流域水资源分配、协调人地关系、调水调沙、防洪减灾等方面都具有重要的理论指导意义。
附件 1 给出了位于小浪底水库下游黄河某水文站近 6 年的水位、水流量与含沙量的实际监测数据，附件 2 给出了该水文站近 6 年黄河断面的测量数据，附件 3 给出了该水文站部分监测点的相关数据。请建立数学模型研究以下问题：
问题 1 研究该水文站黄河水的含沙量与时间、水位、水流量的关系，并估算近 6 年该水文站的年总水流量和年总排沙量。
问题 2 分析近 6 年该水文站水沙通量的突变性、季节性和周期性等特性，研究水沙通量的变化规律。
问题 3 根据该水文站水沙通量的变化规律，预测分析该水文站未来两年水沙通量的变化趋势，并为该水文站制订未来两年最优的采样监测方案（采样监测次数和具体时间等），使其既能及时掌握水沙通量的动态变化情况，又能最大程度地减少监测成本资源。
问题 4 根据该水文站的水沙通量和河底高程的变化情况，分析每年 6-7 月小浪底水库进行“调水调沙”的实际效果。如果不进行“调水调沙”，10 年以后该水文站的河底高程会如何？
附件 1 2016-2021 年黄河水沙监测数据
附件 2 黄河断面的测量数据
附件 3 黄河部分监测点的监测数据
附录 说明
(1) “水位”和“河底高程”均以“1985 国家高程基准”（海拔 72.26 米）为基准面。
(2) 附件中的“起点距离”以河岸边某定点作为起点。
整体求解过程概述(摘要)
黄河是中国的母亲河，其水资源和水沙情况的准确监测对于维护国家生态安全和水资源管理至关重要。本文分析了位于小浪底水库下游黄河某水文站近六年的水位、水流量与含沙量的实际检测数据，基于时间序列预测模型，建立了含沙量预测模型，并分析了该水文站水沙通量的突变性、季节性和周期性等特性，为水文站制定了未来两年最优的采样监测方案，最后分析了“调水调沙”对该水文站河底高程的影响。
针对问题一，本文使用了线性回归模型，来探讨在特定水文站观测到的黄河水的含沙量与时间、水位和水流量之间的关联关系。以描述这些时间、水位和水流量与含沙量之间的关系，并确定各自的系数，从而更好地理解这些因素对含沙量的影响。并给出了水位、流量和含沙量随时间变化的示意图。
针对问题二，采用滑动窗口分析，我们能够识别和量化水沙通量时间序列中的突变点，从而揭示数据中的异常变化。接着，我们进行季节性分解，将时间序列数据分解成长期趋势、季节性成分和周期性成分，以更全面地理解水沙通量的季节性和周期性特征。
针对问题三，根据问题二得出的水沙通量变化规律，使用 ARIMA 模型对未来两年该水文站的水沙通量变化趋势进行预测。采用遗传算法，并结合预测的水沙通量变化趋势，为该水文站制定了未来两年的采样监测方案。
针对问题四，根据 8-5 月（第二年）该水文站的水沙通量和河底高程变化，对6-7 月未进行调水调沙情况进行预测，使用 DID 方法比较其与进行调水调沙之后的差异，以此来分析调水调沙的实际效果，并预测了 10 年后不进行调水调沙情况下河底高程的变化。
模型假设：
1、水位”和“河底高程”均以“1985国家高程基准”（海拔72.26米）为基准面。
2、附件中的“起点距离”以河岸边某定点作为起点。
问题分析：
问题一的分析
问题一要求研究该水文站黄河水的含沙量与时间、水位、水流量的关系，并估算近 6 年该水文站的年总水流量和年总排沙量。本文对附件 1 中该水文站不同时间段下的水位、水流量以及含沙量进行分析，构建含沙量同水位、流量、时间的关系模型，利用已有的含沙量数据对其他时间段含沙量进行预测，根据预测结果来估算该水文站的总水流量和总排沙量。
问题二的分析
问题二要求分析近 6 年该水文站水沙通量的突变性、季节性和周期性等特性，研究水沙通量的变化规律。本文在问题一的基础上，计算每个时间段下的水沙通量，利用滑动窗口分析近六年水沙通量的突变点，使用箱线图直观感受突变点的分布，在对水沙通量进行季节性分解，分析其中的季节性因素和周期性因素，结合上述三点来分析水沙通量的变化规律。
问题三的分析
问题三要求根据该水文站水沙通量的变化规律，预测分析该水文站未来两年水沙通量的变化趋势，并为该水文站制订未来两年最优的采样监测方案，使其既能及时掌握水沙通量的动态变化情况，又能最大程度地减少监测成本资源。本文根据问题二中得到的水沙通量变化规律，使用 ARIMA 模型对水文站未来两年水沙通量变化趋势进行预测。进一步采用遗传算法，结合预测得到的水沙通量变化情况，制定未来两年最优的采样监测方案，通过计算最小成本来尽可能减少监测成本资源。
问题四的分析
问题四要求根据该水文站的水沙通量和河底高程的变化情况，分析每年 6-7 月小浪底水库进行“调水调沙”的实际效果。如果不进行“调水调沙”，10 年以后该水文站的河底高程变化。本文根据问题二中得到的水沙通量，使用月平均采样得到该水文站平均月水沙通量，利用 8 月到次年 5 月的数据对 6-7 月份不进行“调水调沙”情况下的水沙通量进行预测。使用 DID 方法，根据预测得到的水沙通量同实际 6-7 月份水沙通量得到 DID 差异指标，分析“调水调沙”的实际效果。根据附件 2 中同日期下不同起点距离的河底高程，以及附件 3 中的同日期下不同起点距离的水位和水深，计算每日的平均河底高程，进一步计算年平均高程,预测不进行“调水调沙”情况下 10 年后该水文站的河底高程。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：(代码和文档not free)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import pmdarima as pm
from sklearn.linear_model import LinearRegression
imports eaborn as sns
from statsmodels.tsa.arima.model import ARIMA
table=pd.read_excel(r"./data/附件 1.xlsx")
foriinrange(2017,2017+5):
#移除table最后一条数据（重复了）
#print(table.iloc[len(table)-1])
table.drop((len(table)-1),inplace=True)
i=str(i)
temp=pd.read_excel(r"./data/附件 1.xlsx",sheet_name=i)
table=pd.concat([table,temp])
table=table.reset_index(drop=True)
table
#补齐时间
table['年'].fillna(method='ffill',inplace=True)
table['月'].fillna(method='ffill',inplace=True)
table['日'].fillna(method='ffill',inplace=True)
table
#数据预处理
time_list=[]
foriinrange(len(table)):
m,d,h=str(int(table.iloc[i,1])),str(int(table.iloc[i,2])),str(table.iloc[i,3])
if(int(table.iloc[i,1])<10):
m="0"+str(int(table.iloc[i,1]))
if(int(table.iloc[i,2])<10):
d="0"+str(int(table.iloc[i,2]))
#print(m,d)
time=str(int(table.iloc[i,0]))+"-"+m+"-"+d+""+h
#print(time)
time_list.append(time)
temp=pd.DataFrame(time_list,columns=["时刻"])
temp["时刻"]=pd.to_datetime(temp["时刻"])
#temp.to_csv('example3.csv',index=False)
#temp
table1=pd.concat([table,temp],axis=1)
#table
df=table1.iloc[:,[7,4,5,6]]
df.to_csv('example2.csv',index=False)
#将索引转换为日期时间
#df.set_index("时刻",inplace=True)
df
df["时刻"]=pd.to_datetime(df["时刻"])
#将时间序列转换为数值型特征
df1=df.copy()
df1['时刻']=df1['时刻'].apply(lambdax:x.timestamp())
df1
#提取时间、水位、水流量和含沙量的数据
data=df1[pd.notna(df["含沙量(kg/m3)"])]
X=data[['时刻','水位(m)','流量(m3/s)']]
y=data['含沙量(kg/m3)']
y
#建立线性回归模型
model=LinearRegression()
model.fit(X,y)
new_df=df1[pd.isna(df.loc[:,"含沙量(kg/m3)"])]
new_X=new_df.loc[:,['时刻','水位(m)','流量(m3/s)']]
new_df.loc[:,"含沙量(kg/m3)"]=model.predict(new_X)
new_df
#使用fillna方法填充空白部分
table['含沙量(kg/m3)'].fillna(new_df['含沙量(kg/m3)'],inplace=True)
df['含沙量(kg/m3)'].fillna(new_df['含沙量(kg/m3)'],inplace=True)
#table.to_csv('example.csv',index=False)
table
#In[242]:
#计算每年的总水流量和总排沙量
yearly_data=table.groupby(table["年"]).agg({'流量(m3/s)':'sum','含沙量(kg/m3)':'sum'})
#输出近 6 年的年总水流量和年总排沙量
print('近 6 年的年总水流量为：',yearly_data['流量(m3/s)'].sum(),'m³')
print('近 6 年的年总排沙量为：',yearly_data['含沙量(kg/m3)'].sum(),'t')
#In[243]:
#计算水沙通量
df["水沙通量"]=df['含沙量(kg/m3)']*df['流量(m3/s)']
df
#In[14]:
#读取数据
data=pd.read_csv('example2.csv')
#设置日期时间列为索引
data.set_index('时刻',inplace=True)
#创建子图
fig,axes=plt.subplots(nrows=3,ncols=1,figsize=(10,10))
#绘制水位数据
axes[0].plot(data.index,data['水位(m)'],label='WaterLevel',color='blue')
axes[0].set_ylabel('WaterLevel(m)')
axes[0].set_title('WaterLevelOverTime')
#绘制水流量数据
axes[1].plot(data.index,data['流量(m3/s)'],label='FlowRate',color='green')
axes[1].set_ylabel('FlowRate(m^3/s)')
axes[1].set_title('FlowRateOverTime')
#绘制含沙量数据
axes[2].plot(data.index,data['含沙量(kg/m3)'],label='SedimentContent',color='red')
axes[2].set_xlabel('Time')
axes[2].set_ylabel('SedimentContent')
axes[2].set_title('SedimentContentOverTime')
#添加图例
foraxinaxes:
ax.legend()
#调整子图布局
plt.tight_layout()
#显示图形
plt.show()
##分析近 6 年水沙通量的突变性、季节性和周期性等特性
###突变性分析
df
#滑动窗口分析
#定义滑动窗口的大小，这里设置为 10
window_size=10
#创建一个空的DataFrame用于存储突变点
change_points=pd.DataFrame(columns=['时刻','水位(m)','流量(m3/s)','含沙量(kg/m3)','水沙通量'])
#进行滑动窗口分析
foriinrange(len(df)-window_size+1):
window=df.iloc[i:i+window_size]
#计算窗口内数据的均值和标准差
mean_values=window.iloc[:,[4]].mean()
std_values=window.iloc[:,[4]].std()
#设置阈值，可以根据实际情况调整
threshold=2.8#假设阈值为 2
#检测是否有数据超过阈值，如果有，则认为有突变点
if(window.iloc[:,[4]]-mean_values).abs().max().max()>threshold*std_values.max():
cp=pd.DataFrame(window.iloc[-1,:]).T
change_points=pd.concat([change_points,cp])#将突变点添加到结果DataFrame中
#打印突变点
print("突变点:")
print(change_points)
change_points
#创建一个新的Figure
plt.figure(figsize=(12,6))
plt.subplot(411)
plt.boxplot(df['水位(m)'],labels=['waterlevel'],vert=False)
plt.title('waterlevelBoxPlot')
plt.subplot(412)
plt.boxplot(df['流量(m3/s)'],labels=['FlowRate'],vert=False)
plt.title('FlowRateBoxPlot')
plt.subplot(413)
plt.boxplot(df['含沙量(kg/m3)'],labels=['SedimentContent'],vert=False)
plt.title('SedimentContentBoxPlot')
plt.subplot(414)
plt.boxplot(df['水沙通量'],labels=['WaterAndSedimentFlux'],vert=False)
plt.title('WaterAndSedimentFluxBoxPlot')
#显示图形
plt.show()
#创建一个新的Figure
plt.figure(figsize=(12,6))
#可视化水位数据
plt.subplot(311)
plt.plot(df['时刻'],df['水位(m)'],label='waterlevel',color='blue')
plt.xlabel('Time')
plt.ylabel('WaterLevel')
plt.title('WaterLevelOverTime')
#可视化水流量数据
plt.subplot(312)
plt.plot(df['时刻'],df['流量(m3/s)'],label='FlowRate',color='green')
plt.xlabel('Time')
plt.ylabel('FlowRate')
plt.title('FlowRateOverTime')
#可视化含沙量数据
plt.subplot(313)
plt.plot(df['时刻'],df['含沙量(kg/m3)'],label='SedimentContent',color='red')
plt.xlabel('Time')
plt.ylabel('SedimentContent')
plt.title('SedimentContentOverTime')
#在图上标记突变点
forindex,rowinchange_points.iterrows():
plt.subplot(311)
plt.axvline(row['时刻'],color='gray',linestyle='--',linewidth=1)
plt.annotate('change',xy=(row['时刻'],df['水位(m)'].max()),xytext=(-20,30), textcoords='offsetpoints',arrowprops=dict(arrowstyle="->",color='gray'))
plt.subplot(312)
plt.axvline(row['时刻'],color='gray',linestyle='--',linewidth=1)
plt.annotate('change',xy=(row['时刻'],df['流量(m3/s)'].max()),xytext=(-20,30), textcoords='offsetpoints',arrowprops=dict(arrowstyle="->",color='gray'))
lt.subplot(313)
plt.axvline(row['时刻'],color='gray',linestyle='--',linewidth=1)
plt.annotate('change',xy=(row['时刻'],df['含沙量(kg/m3)'].max()),xytext=(-20,30), textcoords='offsetpoints',arrowprops=dict(arrowstyle="->",color='gray'))
#调整子图的布局
plt.tight_layout()
#显示图形
plt.show()
#In[250]:
#将索引转换为日期时间
df.set_index("时刻",inplace=True)
df
#In[255]:
#计算每日季节性成分
seasonal_window=12#每年季节性
seasonal=df.rolling(window=seasonal_window,min_periods=1).mean()
#计算趋势
trend=df-seasonal
#可视化分解结果
plt.figure(figsize=(12,8))
plt.subplot(311)
plt.plot(df['水沙通量'],label='Original')
plt.legend(loc='best')
plt.subplot(312)
plt.plot(trend['水沙通量'],label='Trend')
plt.legend(loc='best')
plt.subplot(313)
plt.plot(seasonal['水沙通量'],label='Seasonal')
plt.legend(loc='best')
plt.tight_layout()
#显示图形
plt.show()