数学建模竞赛中黄河水沙数据的时空特征分析与Python实现

1. 赛题背景与核心任务拆解

每年九月的全国大学生数学建模竞赛，对很多理工科学生来说，都是一场硬仗。2023年的E题“黄河水沙监测数据分析”，直接把战场拉到了黄河水文这个既宏大又具体的领域。题目给的不是抽象的数学公式，而是实打实的黄河干流上四个水文站——小浪底、花园口、高村、利津——从2018年到2022年整整五年的日尺度水沙监测数据。数据列包括日期、水位、流量、含沙量、输沙率，这就是我们手里的全部“弹药”。

这道题的第一问，要求我们“分析水沙通量的时空变化特征及影响因素”。听起来很学术，但翻译成大白话就是：你得从这堆数据里，看出黄河的水和沙是怎么随着时间和空间“动”起来的，并且还得说清楚，为什么它会这么“动”。这不仅仅是画几条折线图、算几个平均值那么简单。它要求你像一个真正的流域分析师一样，从海量、杂乱、可能还有缺失的日数据中，提炼出年际、年内（季节、月份）、甚至日尺度的变化规律，并对比上下游四个站点的异同，最后把这些变化和可能的驱动因素（比如降水、水库调度、人类活动）联系起来。

很多同学一看到“时空变化”、“影响因素”这种词就头疼，觉得无从下手。其实，这道题的魅力恰恰在于它给了你一个非常真实的科研场景：数据是真实的，问题是开放的，没有标准答案，只有逻辑是否自洽、证据是否充分的分析。你的目标不是复现某个经典模型，而是用数据讲一个关于黄河水沙的、令人信服的故事。接下来，我就结合自己处理这类问题的经验，把整个分析链条掰开揉碎，从数据清洗到特征提取，再到可视化与成因推断，一步步带你走完。

2. 数据预处理：从原始表格到规整DataFrame

拿到竞赛数据，第一步永远不是急着跑模型，而是静下心来“盘”数据。原始数据往往以Excel或CSV格式提供，直接pd.read_csv读进来只是万里长征第一步。对于时间序列数据，尤其是水文数据，预处理的质量直接决定了后续所有分析的可靠性。

2.1 数据读取与初步探查

我习惯用pandas配合openpyxl（如果数据是.xlsx格式）来读取数据。读入后，立刻用df.info()和df.head()、df.tail()快速浏览数据结构、数据类型和样本值。

import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 # 假设数据文件为 '黄河水沙数据.csv' df = pd.read_csv('黄河水沙数据.csv', encoding='gbk') # 注意编码，可能是gbk或utf-8 print("数据形状:", df.shape) print("\n数据概览:") print(df.info()) print("\n前5行数据:") print(df.head()) print("\n统计描述:") print(df.describe())

这一步你可能会发现几个典型问题：

1. 列名不规范：原始列名可能是中文，且包含空格或特殊字符。建议统一重命名为英文或简洁中文，方便后续代码引用。例如，将“日期”重命名为date，“流量(m³/s)”重命名为flow。
2. 日期格式混乱：date列可能是字符串格式，如“2018-01-01”、“2018/1/1”或“20180101”。必须用pd.to_datetime()将其统一转换为pandas的datetime类型，这是时间序列分析的基石。
3. 站点信息混合：数据可能是一个大表，包含所有站点的记录，通过“站名”或“站点”列区分。也可能分成了四个独立的文件。我们需要将其整合或拆分，以便进行站点间的对比分析。

2.2 关键清洗步骤：处理缺失值与异常值

水文数据，特别是含沙量和输沙率，受测量条件和极端天气影响，缺失值和异常值非常常见。

处理缺失值：

• 探查：使用df.isnull().sum()统计各列缺失数量。
• 策略：
  • 对于连续缺失少于3天的数据，可以考虑用前后值的线性插值（df.interpolate(method='linear')）填充。这比较符合水文过程的连续性。
  • 对于长时间段（如超过一周）的缺失，或者位于序列开头/结尾的缺失，线性插值可能引入较大误差。这时，一个更稳健的做法是使用该站点同年同月（或相邻日）的历史均值进行填充。例如，计算2018年1月1日所有年份的流量均值，来填充2019年1月1日的缺失值。
  • 切忌：不要简单地用fillna(0)，对于流量、含沙量，零值具有明确的物理意义（断流、无沙），随意填充零会严重扭曲统计结果。

处理异常值：

• 识别：结合统计方法和物理意义。
  • 统计方法：计算每个水文变量的上下四分位数（Q1, Q3）和四分位距（IQR），将大于Q3 + 1.5 * IQR或小于Q1 - 1.5 * IQR的值视为统计异常值。
  • 物理界限：更重要的是物理常识。例如，黄河花园口站的多年平均流量大约在1000-2000 m³/s，如果出现一个50000 m³/s的值，那很可能是数据录入错误或传感器故障，即使它在统计上不算异常，也应视为异常值。
• 处理：
  • 对于明显的录入错误（如小数点错位），可以手动修正。
  • 对于无法解释的极端值，稳妥的做法是将其视为缺失值，然后按上述缺失值处理方法处理。直接删除整行记录要谨慎，可能会破坏时间序列的连续性。

# 示例：处理缺失值和异常值（以流量列为例） def clean_hydrological_data(df, station_name): # 筛选特定站点数据 station_df = df[df['站名'] == station_name].copy() station_df['date'] = pd.to_datetime(station_df['日期']) station_df.set_index('date', inplace=True) station_df.sort_index(inplace=True) # 重命名列 station_df.rename(columns={'流量(m³/s)': 'flow', '含沙量(kg/m³)': 'sediment_concentration', '输沙率(kg/s)': 'sediment_rate'}, inplace=True) # 1. 处理缺失值：线性插值（限制最大连续插值天数为3） for col in ['flow', 'sediment_concentration', 'sediment_rate']: station_df[col] = station_df[col].interpolate(method='linear', limit=3) # 对于插值后仍缺失的，用该月的历史日均值填充（更复杂但更合理） # 这里简化处理，用前向填充 station_df[col].fillna(method='ffill', inplace=True) station_df[col].fillna(method='bfill', inplace=True) # 2. 处理异常值（以流量为例） Q1 = station_df['flow'].quantile(0.25) Q3 = station_df['flow'].quantile(0.75) IQR = Q3 - Q1 lower_bound = Q1 - 1.5 * IQR upper_bound = Q3 + 1.5 * IQR # 将统计异常值设为NaN，然后再次用前后值填充（模拟修正） outlier_mask = (station_df['flow'] < lower_bound) | (station_df['flow'] > upper_bound) station_df.loc[outlier_mask, 'flow'] = np.nan station_df['flow'].interpolate(method='linear', inplace=True) return station_df # 对四个站点分别进行处理 stations = ['小浪底', '花园口', '高村', '利津'] cleaned_data = {} for station in stations: cleaned_data[station] = clean_hydrological_data(df, station)

注意：上面的异常值处理只是一个示例。在实际竞赛中，你需要结合每个站点的历史水文资料，设定更合理的物理阈值。例如，可以参考黄河水利委员会发布的各站点历史极值。

2.3 数据规整与特征工程

清洗后的数据是干净的，但还不够“好用”。为了分析时空变化，我们需要从原始的日数据中，提取出不同时间尺度的特征。

• 时间尺度聚合：

  • 年际变化：计算每年的平均流量、年均含沙量、年总输沙量。这是看大趋势。
  • 年内变化（季节/月度）：计算每个站点每年内各个月份或季节（春季3-5月、夏季6-8月、秋季9-11月、冬季12-2月）的统计量（均值、总量、最大值）。这是看周期性规律。
  • 日尺度波动：可以计算滑动平均（如7日滑动平均）来平滑噪音，观察短期波动趋势。

• 衍生特征计算：

  • 水沙关系：输沙率理论上等于流量乘以含沙量。可以计算sediment_rate_calc = flow * sediment_concentration，并与实测输沙率对比，验证数据的一致性，或发现异常。
  • 峰现时间：找出每年流量和输沙率的峰值（最大值）及其出现的日期，分析“洪峰”和“沙峰”是否同步，这是水沙耦合关系的关键。
  • 变化幅度：计算流量和含沙量的变异系数（标准差/均值），衡量其波动剧烈程度。

# 示例：计算年际和月度特征 def extract_features(station_df, station_name): features = {} df = station_df.copy() # 年际特征 df['year'] = df.index.year annual_flow = df.groupby('year')['flow'].mean() annual_sediment = df.groupby('year')['sediment_rate'].sum() # 年输沙总量 features['annual_flow'] = annual_flow features['annual_sediment'] = annual_sediment # 月度特征（多年平均） df['month'] = df.index.month monthly_flow = df.groupby('month')['flow'].mean() monthly_sediment_concentration = df.groupby('month')['sediment_concentration'].mean() features['monthly_flow'] = monthly_flow features['monthly_sediment_concentration'] = monthly_sediment_concentration # 计算滑动平均（7天） df['flow_7d_avg'] = df['flow'].rolling(window=7, center=True).mean() features['smoothed_flow'] = df['flow_7d_avg'] return features, df # 为每个站点提取特征 features_dict = {} for station in stations: feat, df_with_feat = extract_features(cleaned_data[station], station) features_dict[station] = feat cleaned_data[station] = df_with_feat # 更新数据框，包含新特征

完成这一步，你手里就有了从日到年到站点对比的多维度、规整的数据集，这才是进行深度分析的“弹药库”。

3. 时空变化特征的可视化与分析

有了干净的数据和丰富的特征，接下来就是用图表说话。可视化不是为了好看，而是为了揭示模式、发现异常、支撑论点。针对“时空变化”，我们的可视化策略也要从时间和空间两个维度展开。

3.1 时间变化特征可视化

1. 长期趋势（年际变化）： 绘制四个站点年均流量和年总输沙量的折线图。将四个站点的曲线放在同一张图上，可以直观看出从上游（小浪底）到下游（利津）的变化趋势是否一致。

fig, axes = plt.subplots(2, 1, figsize=(14, 10)) # 年均流量 for station in stations: axes[0].plot(features_dict[station]['annual_flow'].index, features_dict[station]['annual_flow'].values, marker='o', label=station) axes[0].set_title('黄河干流主要水文站年均流量变化 (2018-2022)') axes[0].set_xlabel('年份') axes[0].set_ylabel('流量 (m³/s)') axes[0].legend() axes[0].grid(True, linestyle='--', alpha=0.7) # 年总输沙量 for station in stations: # 注意单位转换，年输沙率总和单位是 kg/s * 秒数，通常转换为万吨 seconds_per_year = 365 * 86400 annual_sediment_10k_tons = features_dict[station]['annual_sediment'] * seconds_per_year / 1e10 axes[1].plot(annual_sediment_10k_tons.index, annual_sediment_10k_tons.values, marker='s', label=station) axes[1].set_title('黄河干流主要水文站年输沙量变化 (2018-2022)') axes[1].set_xlabel('年份') axes[1].set_ylabel('输沙量 (万吨)') axes[1].legend() axes[1].grid(True, linestyle='--', alpha=0.7) plt.tight_layout() plt.show()

从这张图里，你可能发现：2018-2022年间，年均流量是否呈现上升或下降趋势？各站点趋势是否同步？年输沙量的变化趋势与流量是否一致？例如，你可能发现流量变化不大，但输沙量显著减少，这就能引出“水沙关系变化”的讨论点。

2. 周期性规律（年内变化）： 绘制多年平均的月均流量和月均含沙量柱状图或折线图。这能清晰展示黄河水沙的年内分配情况。

fig, axes = plt.subplots(1, 2, figsize=(16, 6)) months = range(1, 13) month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] width = 0.2 x = np.arange(len(months)) for idx, station in enumerate(stations): offset = width * (idx - 1.5) axes[0].bar(x + offset, features_dict[station]['monthly_flow'].values, width, label=station) axes[0].set_title('各站点多年平均月均流量') axes[0].set_xlabel('月份') axes[0].set_ylabel('流量 (m³/s)') axes[0].set_xticks(x) axes[0].set_xticklabels(month_names) axes[0].legend() axes[0].grid(True, axis='y', linestyle='--', alpha=0.7) for idx, station in enumerate(stations): axes[1].plot(months, features_dict[station]['monthly_sediment_concentration'].values, marker='o', label=station) axes[1].set_title('各站点多年平均月均含沙量') axes[1].set_xlabel('月份') axes[1].set_ylabel('含沙量 (kg/m³)') axes[1].set_xticks(months) axes[1].set_xticklabels(month_names) axes[1].legend() axes[1].grid(True, linestyle='--', alpha=0.7) plt.tight_layout() plt.show()

典型的黄河水沙年内分布是：流量和含沙量主要集中在汛期（7-10月），尤其是“七下八上”（七月下旬到八月上旬）的主汛期。从图中可以验证这一规律，并比较不同站点汛期峰值的大小和出现时间的差异。例如，小浪底水库的调节作用可能会使下游花园口的汛期流量过程变得平缓。

3. 短期过程与极端事件： 选取一个典型年份（如2021年），绘制该年四个站点的日流量过程线。可以叠加7日滑动平均线以平滑噪音。

year_to_plot = 2021 fig, ax = plt.subplots(figsize=(15, 8)) for station in stations: df_year = cleaned_data[station][cleaned_data[station].index.year == year_to_plot] ax.plot(df_year.index, df_year['flow'], alpha=0.6, linewidth=1, label=f'{station}日流量') ax.plot(df_year.index, df_year['flow_7d_avg'], linewidth=2, label=f'{station}7日滑动平均') ax.set_title(f'{year_to_plot}年黄河干流主要水文站日流量过程线') ax.set_xlabel('日期') ax.set_ylabel('流量 (m³/s)') ax.legend(loc='upper left', ncol=2) ax.grid(True, linestyle='--', alpha=0.5) # 可以添加汛期阴影背景 ax.axvspan(pd.Timestamp(f'{year_to_plot}-07-01'), pd.Timestamp(f'{year_to_plot}-10-01'), alpha=0.2, color='yellow', label='汛期') plt.xticks(rotation=45) plt.tight_layout() plt.show()

这张图能清晰展示洪水的起涨、峰现、退水过程，以及不同站点洪峰传播的时间差。结合含沙量过程线，可以分析“沙峰”是否滞后于“洪峰”，这是判断河道冲淤状态的重要依据。

3.2 空间变化特征可视化

1. 沿程变化对比： 在同一张图上，用箱型图展示四个站点某个关键指标（如年均流量、汛期平均含沙量）的分布，可以直观看出从上游到下游的空间梯度。

fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 准备数据：各站点所有年份的年均流量列表 annual_flow_list = [] annual_sediment_list = [] for station in stations: annual_flow_list.append(features_dict[station]['annual_flow'].values) # 计算年输沙量列表（单位：万吨） sediment_list = [] for year, total_rate in features_dict[station]['annual_sediment'].items(): sediment_tons = total_rate * 365 * 86400 / 1e10 sediment_list.append(sediment_tons) annual_sediment_list.append(sediment_list) # 绘制箱型图 bp1 = axes[0].boxplot(annual_flow_list, labels=stations, patch_artist=True) axes[0].set_title('各站点年均流量分布对比 (2018-2022)') axes[0].set_ylabel('流量 (m³/s)') axes[0].grid(True, axis='y', linestyle='--', alpha=0.7) bp2 = axes[1].boxplot(annual_sediment_list, labels=stations, patch_artist=True) axes[1].set_title('各站点年输沙量分布对比 (2018-2022)') axes[1].set_ylabel('输沙量 (万吨)') axes[1].grid(True, axis='y', linestyle='--', alpha=0.7) plt.tight_layout() plt.show()

箱型图展示了中位数、四分位数和离散程度。你可能会发现，从小浪底到利津，流量箱型图的中位数可能变化不大（因为黄河下游是地上河，支流汇入少，蒸发渗漏与引水消耗大致平衡），但离散程度（箱子高度）可能增加，说明下游流量受人类调节（如引水灌溉）影响更剧烈。输沙量的中位数则很可能显著下降，反映泥沙沿程淤积。

2. 双变量关系空间对比： 绘制每个站点的流量-含沙量散点图，并放在一起对比。

fig, axes = plt.subplots(2, 2, figsize=(12, 10)) axes = axes.flatten() for idx, station in enumerate(stations): ax = axes[idx] df_station = cleaned_data[station] # 为了图清晰，可以采样，比如每10天取一个点 sample_df = df_station.iloc[::10] scatter = ax.scatter(sample_df['flow'], sample_df['sediment_concentration'], c=sample_df.index.month, cmap='viridis', alpha=0.6, s=10) ax.set_title(f'{station}站 流量-含沙量关系') ax.set_xlabel('流量 (m³/s)') ax.set_ylabel('含沙量 (kg/m³)') ax.grid(True, linestyle='--', alpha=0.5) # 添加颜色条表示月份 plt.colorbar(scatter, ax=ax, label='月份') plt.tight_layout() plt.show()

这张图是分析水沙关系的核心。通常，含沙量随流量增加而增加，但并非简单的线性关系，且存在“滞后环”（涨水段和落水段的点据不重合）。对比四个站点，你可以分析：哪个站点的点据最集中？哪个最分散？点据的分布形态（是否呈带状？是否有明显的拐点？）反映了该河段的水沙输送特性及人类活动（如水库拦沙）的影响。

4. 影响因素分析与综合讨论

可视化揭示了现象，分析部分就要解释原因。题目要求分析“影响因素”，这是一个开放性问题，需要你结合数据特征和地理水文知识进行逻辑推理。以下是一些可能的方向和对应的数据分析方法：

1. 气候变化与降水： 这是最宏观的因素。虽然题目没给降水数据，但你可以从流量过程线的形态进行间接推断。例如，如果某一年所有站点的汛期流量峰值都显著偏低且持续时间短，可以结合公开资料（可简要提及“查阅相关水文年鉴”），推断该年黄河流域可能降水偏少。在论文中，可以定性描述这种关联。

2. 水利工程调节（核心因素）： 小浪底水库是黄河中下游最重要的控制性工程。它的存在会极大地改变下游水沙过程。

• 削峰补枯：对比小浪底站和花园口站的流量过程线。如果花园口的洪峰流量明显低于小浪底，且过程更加平缓，这就是水库削峰的直观证据。可以计算“削峰率”（（入库洪峰-出库洪峰）/入库洪峰）。
• 拦沙作用：计算小浪底站的年均含沙量或输沙量，与下游的花园口、高村、利津站对比。如果小浪底站的含沙量显著高于下游站点，尤其是在非汛期，说明水库拦截了大量泥沙，下泄清水。
• 调水调沙：黄河每年会进行调水调沙。这会在流量过程线上产生非常规的、人为控制的洪峰。你可以观察数据中是否在每年6-7月出现一个流量陡增又陡降的“方波”过程，这很可能就是调水调沙的信号。分析调水调沙期间，下游站点的含沙量如何响应，是同步增加还是滞后增加？

3. 河道冲淤与人类活动：

• 水沙关系变化：分别计算四个站点的流量-含沙量关系式（例如，可以用S = aQ^b进行幂函数拟合，其中S为含沙量，Q为流量）。比较参数a和b的空间差异。下游站点的b值可能变小，意味着流量增加对含沙量增加的贡献率降低，可能反映了河道淤积或河床粗化。
• 输沙量沿程衰减：计算从上游到下游，输沙量的衰减率。例如，（小浪底输沙量 - 利津输沙量）/ 小浪底输沙量。这个衰减率就是泥沙在河道中淤积的比例。可以分析这个比例的年际变化，是趋于稳定还是增大？
• 引水影响：农业灌溉和城市供水会从黄河引水，导致流量沿程减少。可以计算相邻站点间的流量差值（如花园口流量 - 高村流量），观察其季节变化。灌溉季（春灌、冬灌）这个差值可能会变大。

如何在论文中呈现分析： 不要只罗列图表，要用“总-分-总”的结构串联你的发现。

1. 总体概括：首先用一两句话总结2018-2022年黄河水沙通量时空变化的总体特征（如“年均流量相对稳定，但年际波动加大；输沙量呈显著下降趋势，且主要集中于汛期”）。
2. 分点详述：
   • 时间变化：从年际、年内、日尺度分别阐述。例如，“年际尺度上，受2019-2020年降水偏丰影响，各站点流量在2020年出现峰值，随后两年回落”；“年内分配高度不均，约70%的径流量和超过85%的输沙量集中在7-10月汛期”；“日过程显示，小浪底水库调度使下游洪峰过程明显坦化”。
   • 空间变化：从上游到下游对比。例如，“流量沿程变化受引水和蒸发影响，花园口至利津段年均损失约X%；输沙量衰减更为剧烈，小浪底至利津的泥沙淤积比例高达Y%，表明中下游河道仍处于淤积态势”。
   • 水沙关系：分析其时空变异。“上游小浪底站水沙关系相对单一，而下游利津站关系散乱，表明人类活动干扰加剧”。
3. 影响因素综合：将上述特征与可能的原因联系起来。用“数据表明...，这很可能是因为...”的句式。例如，“花园口站汛期含沙量峰值较历史资料降低，且滞后于流量峰值，这与小浪底水库汛期拦沙排浑的调度方式密切相关”。“利津站非汛期出现多次小流量高含沙量事件，可能与河口地区河道疏浚或风暴潮扰动再悬浮有关”。

实操心得：在数学建模论文中，分析部分最忌“两张皮”——图表归图表，文字归文字。一定要做到“图-文-数”三位一体。在描述一个现象时，立即引用对应的图表编号（如“如图3所示”）和关键数据（如“年均流量从2018年的1250 m³/s下降至2022年的980 m³/s”）。你的分析逻辑，就是一根线，把散落的图表和数据珍珠串成一条完整的项链。

5. 代码实现的技巧与避坑指南

虽然题目要求的是数据分析，但一份清晰、可复现、注释良好的Python代码，绝对是论文的加分项。这里分享几个在实现上述分析时容易踩的坑和应对技巧。

1. 时间序列处理的陷阱：

• 时区与格式：pd.to_datetime()非常强大，但遇到“2020-2-30”这种不存在的日期，或者混合了“2020年12月1日”和“2020/12/01”的格式，它会报错或产生NaT。务必先用df['date'].unique()查看日期样本，再用errors='coerce'参数将无法解析的设为NaT，然后检查处理。
• 重采样（Resample）：计算月均值时，不要用groupby('month')，这会混淆不同年份的同一个月。正确做法是先用df.resample('M').mean()重采样为月数据，再进行分析。groupby('month')得到的是所有年份一月份的平均，失去了年际变化信息。
• 滑动窗口：rolling(window=7).mean()计算的是“过去7天（包括当天）的均值”，这会导致序列开头有6个NaN。使用center=True参数可以计算“前后各3.5天”的均值，让滑动平均线与原始数据在时间上对齐，但开头和结尾各有window//2个NaN。

2. 大数据量与计算效率： 五年四个站点的日数据，大约有4站 * 5年 * 365天 ≈ 7300行，不算大。但如果你要进行更复杂的计算（如每个站点逐日的水沙关系拟合），循环嵌套可能会慢。善用pandas的向量化操作和groupby.apply。例如，要计算每个站点每年的流量-含沙量拟合参数，可以：

def fit_power_law(group): # group 是一年的数据 Q = group['flow'].values S = group['sediment_concentration'].values # 移除0值，避免对数无穷大 valid = (Q > 0) & (S > 0) if np.sum(valid) < 10: # 数据点太少不拟合 return pd.Series({'a': np.nan, 'b': np.nan, 'r2': np.nan}) Q_valid = Q[valid] S_valid = S[valid] # 拟合 S = a * Q^b，取对数后线性拟合 ln(S) = ln(a) + b * ln(Q) coeffs = np.polyfit(np.log(Q_valid), np.log(S_valid), 1) b = coeffs[0] ln_a = coeffs[1] a = np.exp(ln_a) # 计算R² S_pred = a * (Q_valid ** b) ss_res = np.sum((S_valid - S_pred) ** 2) ss_tot = np.sum((S_valid - np.mean(S_valid)) ** 2) r2 = 1 - (ss_res / ss_tot) return pd.Series({'a': a, 'b': b, 'r2': r2}) # 对每个站点每年的数据应用拟合 result = df.groupby(['站名', df.index.year]).apply(fit_power_law).reset_index()

3. 可视化美化与信息过载：

• 颜色与样式：使用seaborn的调色板（sns.color_palette()）或matplotlib的tab10、Set2等分类色板，让多系列图表更易区分。线型（-,--,:,-.）和标记（o,s,^,D）也要搭配使用。
• 子图排列：当需要对比多个站点或多个指标时，子图（subplots）比挤在一张图里更清晰。但也要注意子图太多会导致单个图太小。4个站点对比，用2x2的布局通常很合适。
• 图例与标注：图例要清晰，避免遮挡数据。对于关键事件（如洪峰、调水调沙开始日），可以用ax.axvline()添加垂直虚线，并用ax.text()进行标注。
• 避免“图表垃圾”：3D图表、爆炸饼图、过于花哨的背景，在科研图表中都是减分项。保持简洁、清晰、信息密度高才是王道。

4. 结果的可复现性：

• 在代码开头设置随机种子（np.random.seed(42)），虽然本文分析可能用不到随机数，但这是一个好习惯。
• 将数据预处理、特征工程、可视化、分析函数都封装成独立的函数或类，这样代码结构清晰，也便于调试和复用。
• 在关键步骤后，使用df.to_csv('cleaned_data.csv', index=False)保存中间结果。这样如果后续分析出错，不必从头运行耗时的清洗步骤。
• 在Jupyter Notebook中，可以使用%load_ext watermark和%watermark魔法命令记录运行环境的包版本，确保他人能复现你的结果。

最后，记住数学建模竞赛的核心是“建模”和“解决问题”，代码和可视化是工具和语言。你的分析报告（论文）才是最终交付物。代码要服务于清晰的逻辑链条和有力的结论。在论文中，可以精选最有代表性的3-5张核心图表，将其他辅助性图表和完整代码放入附录。整个分析过程，从数据质疑开始，到规律发现，再到机理解释，形成一个闭环，这才是评委们最看重的。