E:\Anaconda\python.exe C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py
生成模拟数据...
开始数据预处理...
数据预处理完成
E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 27969 (\N{CJK UNIFIED IDEOGRAPH-6D41}) missing from font(s) DejaVu Sans.
fig.tight_layout()
E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans.
fig.tight_layout()
E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans.
fig.tight_layout()
E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans.
fig.tight_layout()
E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans.
fig.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 23395 (\N{CJK UNIFIED IDEOGRAPH-5B63}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 33410 (\N{CJK UNIFIED IDEOGRAPH-8282}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 24615 (\N{CJK UNIFIED IDEOGRAPH-6027}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 20998 (\N{CJK UNIFIED IDEOGRAPH-5206}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 35299 (\N{CJK UNIFIED IDEOGRAPH-89E3}) missing from font(s) DejaVu Sans.
plt.tight_layout()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 23395 (\N{CJK UNIFIED IDEOGRAPH-5B63}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 33410 (\N{CJK UNIFIED IDEOGRAPH-8282}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 24615 (\N{CJK UNIFIED IDEOGRAPH-6027}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 20998 (\N{CJK UNIFIED IDEOGRAPH-5206}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 35299 (\N{CJK UNIFIED IDEOGRAPH-89E3}) missing from font(s) DejaVu Sans.
plt.savefig('季节性分解.png', dpi=300)
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 27969 (\N{CJK UNIFIED IDEOGRAPH-6D41}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 23395 (\N{CJK UNIFIED IDEOGRAPH-5B63}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 33410 (\N{CJK UNIFIED IDEOGRAPH-8282}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 24615 (\N{CJK UNIFIED IDEOGRAPH-6027}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 20998 (\N{CJK UNIFIED IDEOGRAPH-5206}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 35299 (\N{CJK UNIFIED IDEOGRAPH-89E3}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans.
plt.show()
C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans.
plt.show()
上面是报错内容
下面是代码
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import pywt
from statsmodels.tsa.seasonal import seasonal_decompose
from scipy import signal
class HydrologicalAnalysis:
def __init__(self, data_path=None):
"""
水文时间序列分析工具
:param data_path: 数据文件路径(CSV格式)
"""
self.data = None
if data_path:
self.load_data(data_path)
def load_data(self, file_path):
"""
加载水文监测数据
:param file_path: CSV文件路径
"""
# 读取数据并处理时间格式
self.data = pd.read_csv(file_path, parse_dates=['时间'], index_col='时间')
# 检查必要字段
required_columns = ['水位(m)', '流量(m3/s)', '含沙量(kg/m3)']
if not all(col in self.data.columns for col in required_columns):
raise ValueError("数据文件缺少必要列:水位(m), 流量(m3/s), 含沙量(kg/m3)")
print(f"成功加载数据:{len(self.data)}条记录")
return self.data
def preprocess_data(self):
"""
数据预处理流程
"""
if self.data is None:
raise ValueError("未加载数据,请先调用load_data方法")
print("开始数据预处理...")
# 1. 缺失值处理(前向填充+线性插值)
self.data = self.data.ffill().interpolate(method='linear')
# 2. 异常值处理(3σ原则)
for col in ['流量(m3/s)', '含沙量(kg/m3)']:
mean = self.data[col].mean()
std = self.data[col].std()
self.data[col] = np.where(
(self.data[col] < mean - 3 * std) | (self.data[col] > mean + 3 * std),
mean,
self.data[col]
)
# 3. 重采样为日数据
daily_data = self.data.resample('D').mean()
print("数据预处理完成")
return daily_data
def seasonal_decomposition(self, series, period=12):
"""
时间序列季节性分解
:param series: 时间序列数据
:param period: 季节性周期(月数据默认为12)
:return: 分解结果
"""
result = seasonal_decompose(series, model='additive', period=period)
return result
def wavelet_analysis(self, series, title='水文序列'):
"""
小波分析主函数
:param series: 时间序列数据
:param title: 分析标题
:return: (小波系数, 主周期)
"""
# 1. 参数设置
scales = np.arange(1, 365) # 1天到1年尺度
wavelet = 'morl' # Morlet小波
# 2. 小波变换
coef, freqs = pywt.cwt(series, scales, wavelet)
# 3. 小波方差计算
variance = np.mean(np.abs(coef) ** 2, axis=1)
main_scale = scales[np.argmax(variance)]
# 4. 可视化
self.plot_wavelet_results(series, scales, coef, variance, main_scale, title)
return coef, main_scale
def plot_wavelet_results(self, series, scales, coef, variance, main_scale, title):
"""
绘制小波分析结果
"""
plt.figure(figsize=(15, 12))
# 1. 原始序列图
plt.subplot(3, 1, 1)
plt.plot(series.index, series)
plt.title(f'{title} - 原始时间序列')
plt.xlabel('日期')
plt.ylabel('数值')
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
# 2. 小波系数实部等值线图
plt.subplot(3, 1, 2)
plt.contourf(series.index, scales, np.real(coef), cmap='jet', levels=100)
plt.colorbar(label='小波系数实部')
plt.title(f'{title} - 小波系数实部等值线图')
plt.ylabel('时间尺度(天)')
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
# 3. 小波方差图
plt.subplot(3, 1, 3)
plt.plot(scales, variance)
plt.axvline(main_scale, color='red', linestyle='--',
label=f'主周期: {main_scale}天')
plt.title(f'{title} - 小波方差分析')
plt.xlabel('时间尺度(天)')
plt.ylabel('方差')
plt.legend()
plt.tight_layout()
plt.savefig(f'{title}_小波分析.png', dpi=300)
plt.show()
print(f"{title}主周期: {main_scale}天")
def full_analysis(self):
"""
完整分析流程
"""
# 1. 数据预处理
daily_data = self.preprocess_data()
# 2. 季节性分解
plt.figure(figsize=(15, 10))
for i, col in enumerate(['流量(m3/s)', '含沙量(kg/m3)'], 1):
plt.subplot(2, 1, i)
result = self.seasonal_decomposition(daily_data[col])
result.plot()
plt.title(f'{col}季节性分解')
plt.tight_layout()
plt.savefig('季节性分解.png', dpi=300)
plt.show()
# 3. 小波分析
flow_coef, flow_period = self.wavelet_analysis(
daily_data['流量(m3/s)'], '流量'
)
sediment_coef, sediment_period = self.wavelet_analysis(
daily_data['含沙量(kg/m3)'], '含沙量'
)
# 4. 交叉小波分析(流量与含沙量关系)
self.cross_wavelet_analysis(
daily_data['流量(m3/s)'],
daily_data['含沙量(kg/m3)'],
'流量-含沙量'
)
return {
'flow_period': flow_period,
'sediment_period': sediment_period
}
def cross_wavelet_analysis(self, series1, series2, title='交叉分析'):
"""
交叉小波分析
:param series1: 第一个时间序列
:param series2: 第二个时间序列
:param title: 分析标题
"""
# 1. 计算小波变换
scales = np.arange(1, 365)
wavelet = 'morl'
coef1, _ = pywt.cwt(series1, scales, wavelet)
coef2, _ = pywt.cwt(series2, scales, wavelet)
# 2. 计算交叉小波谱
cross_spectrum = coef1 * np.conj(coef2)
# 3. 可视化
plt.figure(figsize=(15, 8))
# 交叉小波谱实部
plt.subplot(2, 1, 1)
plt.contourf(series1.index, scales, np.real(cross_spectrum),
cmap='RdBu_r', levels=100)
plt.colorbar(label='实部')
plt.title(f'{title} - 交叉小波谱实部')
plt.ylabel('时间尺度(天)')
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
# 交叉小波谱相位
plt.subplot(2, 1, 2)
phase = np.angle(cross_spectrum)
plt.contourf(series1.index, scales, phase, cmap='hsv', levels=100)
plt.colorbar(label='相位(弧度)')
plt.title(f'{title} - 交叉小波谱相位')
plt.xlabel('日期')
plt.ylabel('时间尺度(天)')
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
plt.tight_layout()
plt.savefig(f'{title}_交叉小波分析.png', dpi=300)
plt.show()
# ======================
# 使用示例
# ======================
if __name__ == "__main__":
# 1. 创建分析对象
analyzer = HydrologicalAnalysis()
# 2. 加载数据(替换为实际文件路径)
# analyzer.load_data('水文监测数据.csv')
# 3. 生成模拟数据(实际使用时请注释掉)
print("生成模拟数据...")
dates = pd.date_range(start='2016-01-01', end='2021-12-31', freq='D')
flow = np.sin(2 * np.pi * dates.dayofyear / 365) * 100 + 500 + np.random.normal(0, 50, len(dates))
sediment = np.cos(2 * np.pi * dates.dayofyear / 365) * 2 + 5 + np.random.normal(0, 1, len(dates))
analyzer.data = pd.DataFrame({
'水位(m)': np.random.uniform(40, 45, len(dates)),
'流量(m3/s)': flow,
'含沙量(kg/m3)': sediment
}, index=dates)
# 4. 执行完整分析
results = analyzer.full_analysis()
print("\n分析结果摘要:")
print(f"流量主周期: {results['flow_period']}天")
print(f"含沙量主周期: {results['sediment_period']}天")