import matplotlib.pyplot as plt import numpy as np import matplotlib import pandas as pd from datetime import datetime, timedelta from scipy.signal import savgol_filter # 设置中文字体支持 matplotlib.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 matplotlib.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 # 定义改进的盐度数据异常过滤方法 def filter_salinity_anomalies(df, threshold_ratio=0.5, window_size=5, max_days=1): # 复制数据，避免修改原始数据 filtered_df = df.copy() values = filtered_df['Value'].values dates = filtered_df['DateTime'].values # 1. 首先处理单个异常点 i = 1 while i < len(values): # 检查当前值是否小于前一个值的threshold_ratio if values[i] < values[i-1] * threshold_ratio: baseline = values[i-1] # 基准值为上一个正常的盐度值 anomaly_start = i j = i # 向后查找，直到找到一个不小于基准值threshold_ratio的点 # 或者直到时间区间超过max_days天 anomaly_start_date = dates[anomaly_start] max_date = anomaly_start_date + np.timedelta64(int(max_days*24), 'h') while j < len(values) and values[j] < baseline * threshold_ratio and dates[j] <= max_date: j += 1 anomaly_end = j - 1 # 异常区间的结束位置 # 处理异常区间 if anomaly_end - anomaly_start < 3: # 短区间用线性插值 if j < len(values): # 如果异常区间后还有数据点，使用线性插值 for k in range(anomaly_start, anomaly_end + 1): # 线性插值：在基准值和异常区间后第一个正常值之间进行平滑过渡 ratio = (k - anomaly_start + 1) / (anomaly_end - anomaly_start + 2) values[k] = baseline * (1 - ratio) + values[j] * ratio # 确保平滑后的值不低于基准的threshold_ratio values[k] = max(values[k], baseline * threshold_ratio) else: # 如果异常区间到数据末尾，使用基准值的threshold_ratio填充 for k in range(anomaly_start, anomaly_end + 1): values[k] = baseline * threshold_ratio else: # 长区间使用更简单的平滑方式，避免插值错误 # 使用线性插值来避免非有限值问题 if j < len(values): end_val = values[j] # 为每个点创建线性插值 for k in range(anomaly_start, anomaly_end + 1): fraction = (k - anomaly_start) / (j - anomaly_start) if j > anomaly_start else 0 interpolated = baseline * (1 - fraction) + end_val * fraction values[k] = max(interpolated, baseline * threshold_ratio) else: # 如果异常区间到数据末尾，使用基准值的threshold_ratio填充 for k in range(anomaly_start, anomaly_end + 1): values[k] = baseline * threshold_ratio i = j # 跳过已处理的异常区间 else: i += 1 # 2. 应用Savitzky-Golay滤波进行整体平滑 if len(values) > window_size: # 确保window_size是奇数 if window_size % 2 == 0: window_size += 1 # 应用Savitzky-Golay滤波 try: # 对数据进行平滑，但保留原始的特性 smoothed = savgol_filter(values, window_size, 3) # 确保平滑后的数据不会小于相邻点的threshold_ratio for i in range(1, len(smoothed)): smoothed[i] = max(smoothed[i], smoothed[i-1] * threshold_ratio) values = smoothed except Exception as e: print(f"Savitzky-Golay滤波应用失败: {e}") filtered_df['Value'] = values return filtered_df # 读取CSV数据 df1 = pd.read_csv('yuce_data/一取水.csv') df2 = pd.read_csv('yuce_data/青龙港1.csv') df3 = pd.read_csv('yuce_data/太仓石化盐度2.csv') # 将时间列转换为datetime格式 df1['DateTime'] = pd.to_datetime(df1['DateTime']) df2['DateTime'] = pd.to_datetime(df2['DateTime']) df3['DateTime'] = pd.to_datetime(df3['DateTime']) # 应用盐度数据异常过滤方法 filtered_df1 = filter_salinity_anomalies(df1, threshold_ratio=0.5, window_size=7, max_days=1) filtered_df2 = filter_salinity_anomalies(df2, threshold_ratio=0.5, window_size=7, max_days=1) filtered_df3 = filter_salinity_anomalies(df3, threshold_ratio=0.5, window_size=7, max_days=1) # 创建图表，调整为4个子图 plt.figure(figsize=(16, 16)) # 绘制一取水的原始数据和过滤后的数据 plt.subplot(4, 1, 1) plt.plot(df1['DateTime'], df1['Value'], 'o-', color='lightblue', linewidth=1, markersize=2, alpha=0.6, label='一取水(原始)') plt.plot(filtered_df1['DateTime'], filtered_df1['Value'], 'o-', color='blue', linewidth=1.5, markersize=3, label='一取水(过滤后)') plt.title('一取水盐度监测数据异常过滤', fontsize=16) plt.ylabel('盐度值', fontsize=14) plt.legend(loc='best', fontsize=10) plt.grid(True, linestyle='--', alpha=0.7) # 绘制青龙港1的原始数据和过滤后的数据 plt.subplot(4, 1, 2) plt.plot(df2['DateTime'], df2['Value'], 's-', color='lightcoral', linewidth=1, markersize=2, alpha=0.6, label='青龙港1(原始)') plt.plot(filtered_df2['DateTime'], filtered_df2['Value'], 's-', color='red', linewidth=1.5, markersize=3, label='青龙港1(过滤后)') plt.title('青龙港1盐度监测数据异常过滤', fontsize=16) plt.ylabel('盐度值', fontsize=14) plt.legend(loc='best', fontsize=10) plt.grid(True, linestyle='--', alpha=0.7) # 绘制太仓石化盐度1的原始数据和过滤后的数据 plt.subplot(4, 1, 3) plt.plot(df3['DateTime'], df3['Value'], '^-', color='lightgreen', linewidth=1, markersize=2, alpha=0.6, label='太仓石化1(原始)') plt.plot(filtered_df3['DateTime'], filtered_df3['Value'], '^-', color='green', linewidth=1.5, markersize=3, label='太仓石化1(过滤后)') plt.title('太仓石化盐度1监测数据异常过滤', fontsize=16) plt.ylabel('盐度值', fontsize=14) plt.legend(loc='best', fontsize=10) plt.grid(True, linestyle='--', alpha=0.7) # 计算过滤前后的差异 df1['Difference'] = filtered_df1['Value'] - df1['Value'] df2['Difference'] = filtered_df2['Value'] - df2['Value'] df3['Difference'] = filtered_df3['Value'] - df3['Value'] # 绘制差异 plt.subplot(4, 1, 4) plt.plot(df1['DateTime'], df1['Difference'], 'o-', color='blue', linewidth=1.5, markersize=3, label='一取水(差异)') plt.plot(df2['DateTime'], df2['Difference'], 's-', color='red', linewidth=1.5, markersize=3, label='青龙港1(差异)') plt.plot(df3['DateTime'], df3['Difference'], '^-', color='green', linewidth=1.5, markersize=3, label='太仓石化1(差异)') plt.title('数据处理前后差异', fontsize=16) plt.xlabel('时间', fontsize=14) plt.ylabel('差异值', fontsize=14) plt.legend(loc='best', fontsize=10) plt.grid(True, linestyle='--', alpha=0.7) # 格式化x轴日期显示 for ax in plt.gcf().get_axes(): ax.xaxis_date() plt.gcf().autofmt_xdate() # 显示图表 plt.tight_layout() plt.show() # 保存过滤后的数据 filtered_df1.to_csv('yuce_data/一取水_filtered.csv', index=False) filtered_df2.to_csv('yuce_data/青龙港1_filtered.csv', index=False) filtered_df3.to_csv('yuce_data/太仓石化盐度1_filtered.csv', index=False) # 打印过滤前后的统计信息及处理区间 print("\n数据过滤前后统计信息:") print("\n一取水:") print(f"原始数据 - 最小值: {df1['Value'].min():.2f}, 最大值: {df1['Value'].max():.2f}, 均值: {df1['Value'].mean():.2f}, 标准差: {df1['Value'].std():.2f}") print(f"过滤后数据 - 最小值: {filtered_df1['Value'].min():.2f}, 最大值: {filtered_df1['Value'].max():.2f}, 均值: {filtered_df1['Value'].mean():.2f}, 标准差: {filtered_df1['Value'].std():.2f}") print("\n青龙港1:") print(f"原始数据 - 最小值: {df2['Value'].min():.2f}, 最大值: {df2['Value'].max():.2f}, 均值: {df2['Value'].mean():.2f}, 标准差: {df2['Value'].std():.2f}") print(f"过滤后数据 - 最小值: {filtered_df2['Value'].min():.2f}, 最大值: {filtered_df2['Value'].max():.2f}, 均值: {filtered_df2['Value'].mean():.2f}, 标准差: {filtered_df2['Value'].std():.2f}") print("\n太仓石化盐度1:") print(f"原始数据 - 最小值: {df3['Value'].min():.2f}, 最大值: {df3['Value'].max():.2f}, 均值: {df3['Value'].mean():.2f}, 标准差: {df3['Value'].std():.2f}") print(f"过滤后数据 - 最小值: {filtered_df3['Value'].min():.2f}, 最大值: {filtered_df3['Value'].max():.2f}, 均值: {filtered_df3['Value'].mean():.2f}, 标准差: {filtered_df3['Value'].std():.2f}") # 找出异常点，计算处理的区间数量和时间跨度 def get_anomaly_stats(df_orig, df_filtered): diff = abs(df_filtered['Value'] - df_orig['Value']) anomalies = df_orig[diff > 0.01] # 差异大于0.01的点认为是异常点 if len(anomalies) == 0: return 0, "无异常区间" # 计算连续的异常区间 anomaly_indices = anomalies.index.tolist() intervals = [] start_idx = anomaly_indices[0] prev_idx = start_idx for idx in anomaly_indices[1:]: if idx > prev_idx + 1: # 不连续 intervals.append((start_idx, prev_idx)) start_idx = idx prev_idx = idx intervals.append((start_idx, prev_idx)) # 添加最后一个区间 # 计算区间的时间跨度 time_spans = [] for start, end in intervals: if start == end: span = "单点" else: delta = df_orig.iloc[end]['DateTime'] - df_orig.iloc[start]['DateTime'] hours = delta.total_seconds() / 3600 if hours < 1: span = f"{int(hours * 60)}分钟" elif hours < 24: span = f"{int(hours)}小时" else: span = f"{int(hours/24)}天{int(hours%24)}小时" time_spans.append(span) return len(intervals), time_spans anomaly_count1, time_spans1 = get_anomaly_stats(df1, filtered_df1) anomaly_count2, time_spans2 = get_anomaly_stats(df2, filtered_df2) anomaly_count3, time_spans3 = get_anomaly_stats(df3, filtered_df3) print("\n异常区间统计:") print(f"一取水 - 异常区间数量: {anomaly_count1}") if anomaly_count1 > 0: print(f"区间时间跨度: {time_spans1}") print(f"青龙港1 - 异常区间数量: {anomaly_count2}") if anomaly_count2 > 0: print(f"区间时间跨度: {time_spans2}") print(f"太仓石化盐度1 - 异常区间数量: {anomaly_count3}") if anomaly_count3 > 0: print(f"区间时间跨度: {time_spans3}")