我正在做一个桥梁温致挠度去除与预测的算法，为了在实际应用中表现优异，请参照以下代码写一个将EMD改为VMD的算法，要求选用适合VMD参数自适应的方法，并考虑计算速度import os
import json
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD
from sklearn.ensemble import IsolationForest
from scipy.signal import savgol_filter, find_peaks
from scipy.fft import fft, fftfreq
from sklearn.metrics import mean_squared_error
import warnings

忽略不必要的警告
warnings.filterwarnings("ignore")

解决中文显示问题
plt.rcParams['font.sans-serif'] = ['SimHei'] # 使用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题

参数文件路径
PARAM_FILE = 'best_params.json'

def calculate_dominant_frequency(imf, sampling_rate):
N = len(imf)
yf = fft(imf)
xf = fftfreq(N, d=1/sampling_rate)
positive_freqs = xf[:N//2]
positive_yf = np.abs(yf[:N//2])


if len(positive_yf) == 0:
    return 0.0

idx = np.argmax(positive_yf)
dominant_freq = positive_freqs[idx]
return dominant_freq
def emd_spectral_decomposition(signal, time_seconds, frequency_threshold=None):
emd = EMD()
IMFs = emd.emd(S=signal, T=time_seconds)
num_IMFs = IMFs.shape[0]
print(f"分解得到的IMF数量: {num_IMFs}")


time_deltas = np.diff(time_seconds)
avg_time_delta = np.mean(time_deltas)
if avg_time_delta == 0:
    raise ValueError("时间序列的平均时间间隔为0，无法计算采样率。")
sampling_rate = 1 / avg_time_delta
print(f"估计采样率: {sampling_rate:.4f} Hz")

dominant_freqs = []
for i, imf in enumerate(IMFs):
    dom_freq = calculate_dominant_frequency(imf, sampling_rate)
    dominant_freqs.append(dom_freq)
    print(f"IMF {i+1}: 主频率 = {dom_freq:.4f} Hz")

if frequency_threshold is None:
    frequency_threshold = adaptive_frequency_threshold(dominant_freqs)
    print(f"自动确定的频率阈值为: {frequency_threshold:.4f} Hz")

low_freq_IMF_indices = [i for i, freq in enumerate(dominant_freqs) if freq < frequency_threshold]
print(f"需要去除的低频IMF编号: {low_freq_IMF_indices}")

if low_freq_IMF_indices:
    removed_signal = np.sum(IMFs[low_freq_IMF_indices, :], axis=0)
    cleaned_signal = signal - removed_signal
else:
    cleaned_signal = signal

return cleaned_signal, low_freq_IMF_indices, IMFs, frequency_threshold
def adaptive_frequency_threshold(dominant_freqs):
hist, bin_edges = np.histogram(dominant_freqs, bins='auto', density=True)
threshold_idx = np.argmax(hist)
frequency_threshold = (bin_edges[threshold_idx] + bin_edges[threshold_idx + 1]) / 2
return frequency_threshold

def dynamic_window_size(cleaned_signal):
autocorrelation = np.correlate(cleaned_signal, cleaned_signal, mode='full')
autocorrelation = autocorrelation[autocorrelation.size // 2:]


half_max = autocorrelation[0] * 0.5
peaks, _ = find_peaks(-autocorrelation, height=-half_max)
if len(peaks) > 0:
    window_size = peaks[0]
else:
    window_size = min(24, len(cleaned_signal) // 2)
return max(12, window_size)
def moving_average_forecast(series, window_size):
forecast = []
for t in range(len(series) - window_size + 1):
forecast.append(series[t:t + window_size].mean())
return np.array(forecast)

def plot_imfs(time, IMFs, low_freq_indices):
num_IMFs = IMFs.shape[0]
plt.figure(figsize=(12, 3*num_IMFs))


for i in range(num_IMFs):
    plt.subplot(num_IMFs, 1, i + 1)
    plt.plot(time, IMFs[i], label=f'IMF {i+1}', color='b' if i not in low_freq_indices else 'r')
    if i in low_freq_indices:
        plt.legend(loc='upper right')
        plt.title(f'IMF {i+1} (低频, 温度影响)')
    else:
        plt.legend(loc='upper right')
        plt.title(f'IMF {i+1}')
    plt.xlabel('时间')
    plt.ylabel('挠度')

plt.tight_layout()
plt.show()
def save_best_params(params):
with open(PARAM_FILE, 'w', encoding='utf-8') as f:
json.dump(params, f, ensure_ascii=False, indent=4)
print(f"最佳参数已保存到 {PARAM_FILE}。")

def load_best_params():
if not os.path.exists(PARAM_FILE):
print(f"参数文件 {PARAM_FILE} 不存在。")
return None
with open(PARAM_FILE, 'r', encoding='utf-8') as f:
params = json.load(f)
print(f"已加载最佳参数: {params}")
return params

def estimate_contamination(series):
mean = np.mean(series)
std = np.std(series)
threshold_upper = mean + 3 * std
threshold_lower = mean - 3 * std
num_anomalies = np.sum((series > threshold_upper) | (series < threshold_lower))
contamination = num_anomalies / len(series)
contamination = min(max(contamination, 0.01), 0.05)
print(f"自动估计的异常比例（contamination）: {contamination:.4f}")
return contamination

def main_automated_process(file_path):
try:
df = pd.read_excel(file_path)
except Exception as e:
print(f"读取Excel文件时出错: {e}")
return None, None


print("Excel文件读取成功。")
print(f"数据总行数: {len(df)}")
print(f"列名: {df.columns.tolist()}")

required_columns = ['时间', '挠度']
for col in required_columns:
    if col not in df.columns:
        print(f"Excel文件中缺少 '{col}' 列，请检查列名。")
        return None, None

try:
    df['时间'] = pd.to_datetime(df['时间'])
except Exception as e:
    print(f"转换 '时间' 列为 datetime 格式时出错: {e}")
    return None, None

window_size = 11
poly_order = 3
if window_size % 2 == 0:
    window_size += 1
    print(f"window_size为偶数，已调整为 {window_size}")

try:
    df['平滑挠度'] = savgol_filter(df['挠度'], window_size, poly_order)
except Exception as e:
    print(f"Savgol 滤波器应用时出错: {e}")
    return None, None

contamination = estimate_contamination(df['平滑挠度'].values)

iso_forest = IsolationForest(contamination=contamination, random_state=42)
try:
    anomalies = iso_forest.fit_predict(df[['平滑挠度']])
except Exception as e:
    print(f"孤立森林异常检测时出错: {e}")
    return None, None

df['挠度_去除异常'] = df['平滑挠度']
df.loc[anomalies == -1, '挠度_去除异常'] = np.nan
num_anomalies = (anomalies == -1).sum()
print(f"检测到的异常值数量: {num_anomalies}")

cleaned_data = df[['时间', '挠度_去除异常']].dropna()

if cleaned_data.empty:
    print("所有数据点都被标记为异常，无法继续处理。")
    return df, None

try:
    time_seconds = (cleaned_data['时间'] - cleaned_data['时间'].iloc[0]).dt.total_seconds().values
except Exception as e:
    print(f"转换时间列为秒数时出错: {e}")
    return None, None

try:
    cleaned_signal, removed_IMFs, IMFs, frequency_threshold = emd_spectral_decomposition(
        signal=cleaned_data['挠度_去除异常'].values,
        time_seconds=time_seconds
    )
except Exception as e:
    print(f"EMD分解过程中出错: {e}")
    return df, None

plot_imfs(cleaned_data['时间'], IMFs, removed_IMFs)

df.loc[cleaned_data.index, '去除温致挠度'] = cleaned_signal

window_size_forecast = dynamic_window_size(cleaned_signal)
print(f"动态确定的移动平均窗口大小 (window_size_forecast): {window_size_forecast}")

try:
    forecast = moving_average_forecast(cleaned_signal, window_size_forecast)
    if len(forecast) == 0:
        print("移动平均预测生成的预测数据为空。")
    else:
        actual = cleaned_signal[window_size_forecast -1 :]
        predicted = forecast
        rmse = np.sqrt(mean_squared_error(actual, predicted))
        print(f"移动平均预测的RMSE: {rmse:.4f}")
        forecast_series = np.concatenate([cleaned_signal[:len(cleaned_signal) - len(forecast)], forecast])
        df.loc[cleaned_data.index, '预测'] = forecast_series
except Exception as e:
    print(f"移动平均预测过程中出错: {e}")

used_params = {
    'contamination': contamination,
    'frequency_threshold': frequency_threshold,
    'window_size_forecast': window_size_forecast
}
save_best_params(used_params)

try:
    plt.figure(figsize=(14, 8))
    
    plt.subplot(2, 1, 1)
    plt.plot(cleaned_data['时间'], cleaned_data['挠度_去除异常'], label='原始信号', color='b', alpha=0.5)
    plt.plot(cleaned_data['时间'], cleaned_signal, label='去除低频IMF后的信号', color='g')
    plt.legend(loc='upper right')
    plt.title('原始信号 vs 去除低频IMF后的信号')
    plt.xlabel('时间')
    plt.ylabel('挠度')
    
    if forecast_series.size > 0:
        plt.subplot(2, 1, 2)
        plt.plot(cleaned_data['时间'], cleaned_signal, label='去除低频IMF后的信号', color='g')
        plt.plot(cleaned_data['时间'], forecast_series, label='移动平均预测', color='r')
        plt.legend(loc='upper right')
        plt.title('去除低频IMF后的信号及其预测')
        plt.xlabel('时间')
        plt.ylabel('挠度')
    else:
        plt.subplot(2, 1, 2)
        plt.plot(cleaned_data['时间'], cleaned_signal, label='去除低频IMF后的信号', color='g')
        plt.legend(loc='upper right')
        plt.title('去除低频IMF后的信号')
        plt.xlabel('时间')
        plt.ylabel('挠度')
    
    plt.tight_layout()
    plt.show()
except Exception as e:
    print(f"可视化过程中出错: {e}")

return df, used_params
def load_and_process_with_best_params(file_path, best_params):
contamination = best_params.get('contamination', 0.02)
frequency_threshold = best_params.get('frequency_threshold', 0.1)
window_size_forecast = best_params.get('window_size_forecast', 24)


try:
    df = pd.read_excel(file_path)
except Exception as e:
    print(f"读取Excel文件时出错: {e}")
    return None

print("Excel文件读取成功。")
print(f"数据总行数: {len(df)}")
print(f"列名: {df.columns.tolist()}")

required_columns = ['时间', '挠度']
for col in required_columns:
    if col not in df.columns:
        print(f"Excel文件中缺少 '{col}' 列，请检查列名。")
        return None

try:
    df['时间'] = pd.to_datetime(df['时间'])
except Exception as e:
    print(f"转换 '时间' 列为 datetime 格式时出错: {e}")
    return None

window_size = 11
poly_order = 3
if window_size % 2 == 0:
    window_size += 1
    print(f"window_size为偶数，已调整为 {window_size}")

try:
    df['平滑挠度'] = savgol_filter(df['挠度'], window_size, poly_order)
except Exception as e:
    print(f"Savgol 滤波器应用时出错: {e}")
    return None

iso_forest = IsolationForest(contamination=contamination, random_state=42)
try:
    anomalies = iso_forest.fit_predict(df[['平滑挠度']])
except Exception as e:
    print(f"孤立森林异常检测时出错: {e}")
    return None

df['挠度_去除异常'] = df['平滑挠度']
df.loc[anomalies == -1, '挠度_去除异常'] = np.nan
num_anomalies = (anomalies == -1).sum()
print(f"检测到的异常值数量: {num_anomalies}")

cleaned_data = df[['时间', '挠度_去除异常']].dropna()

if cleaned_data.empty:
    print("所有数据点都被标记为异常，无法继续处理。")
    return df

try:
    time_seconds = (cleaned_data['时间'] - cleaned_data['时间'].iloc[0]).dt.total_seconds().values
except Exception as e:
    print(f"转换时间列为秒数时出错: {e}")
    return None

try:
    cleaned_signal, removed_IMFs, IMFs, _ = emd_spectral_decomposition(
        signal=cleaned_data['挠度_去除异常'].values,
        time_seconds=time_seconds,
        frequency_threshold=frequency_threshold
    )
except Exception as e:
    print(f"EMD分解过程中出错: {e}")
    return df

plot_imfs(cleaned_data['时间'], IMFs, removed_IMFs)

df.loc[cleaned_data.index, '去除温致挠度'] = cleaned_signal

try:
    forecast = moving_average_forecast(cleaned_signal, window_size_forecast)
    if len(forecast) == 0:
        print("移动平均预测生成的预测数据为空。")
    else:
        actual = cleaned_signal[window_size_forecast -1 :]
        predicted = forecast
        rmse = np.sqrt(mean_squared_error(actual, predicted))
        print(f"移动平均预测的RMSE: {rmse:.4f}")
        forecast_series = np.concatenate([cleaned_signal[:len(cleaned_signal) - len(forecast)], forecast])
        df.loc[cleaned_data.index, '预测'] = forecast_series
except Exception as e:
    print(f"移动平均预测过程中出错: {e}")

try:
    plt.figure(figsize=(14, 8))
    
    plt.subplot(2, 1, 1)
    plt.plot(cleaned_data['时间'], cleaned_data['挠度_去除异常'], label='原始信号', color='b', alpha=0.5)
    plt.plot(cleaned_data['时间'], cleaned_signal, label='去除低频IMF后的信号', color='g')
    plt.legend(loc='upper right')
    plt.title('原始信号 vs 去除低频IMF后的信号')
    plt.xlabel('时间')
    plt.ylabel('挠度')
    
    if forecast_series.size > 0:
        plt.subplot(2, 1, 2)
        plt.plot(cleaned_data['时间'], cleaned_signal, label='去除低频IMF后的信号', color='g')
        plt.plot(cleaned_data['时间'], forecast_series, label='移动平均预测', color='r')
        plt.legend(loc='upper right')
        plt.title('去除低频IMF后的信号及其预测')
        plt.xlabel('时间')
        plt.ylabel('挠度')
    else:
        plt.subplot(2, 1, 2)
        plt.plot(cleaned_data['时间'], cleaned_signal, label='去除低频IMF后的信号', color='g')
        plt.legend(loc='upper right')
        plt.title('去除低频IMF后的信号')
        plt.xlabel('时间')
        plt.ylabel('挠度')
    
    plt.tight_layout()
    plt.show()
except Exception as e:
    print(f"可视化过程中出错: {e}")

return df, best_params
def display_menu():
print("\n请选择操作模式:")
print("1. 自适应选择参数 (自动参数优化)")
print("2. 使用上次的参数")
print("3. 退出")
choice = input("请输入选项 (1/2/3): ").strip()
return choice

def main():
file_path = './挠度数据910.xlsx' # 替换为实际文件路径


while True:
    choice = display_menu()
    
    if choice == '1':
        print("\n您选择了自适应选择参数 (自动参数优化)。")
        df_result, optimal_params = main_automated_process(file_path)
        if optimal_params:
            print("\n自动化处理完成。")
            print(f"使用的参数组合: {optimal_params}")
    elif choice == '2':
        print("\n您选择了使用上次的参数。")
        best_params = load_best_params()
        if best_params is None:
            print("未找到上次的参数，请先进行自适应选择参数 (选项1)。")
            continue
        df_result, _ = load_and_process_with_best_params(file_path, best_params)
        if df_result is not None:
            print("\n使用上次的参数处理完成。")
    elif choice == '3':
        print("退出程序。")
        break
    else:
        print("无效的选项，请重新选择。")
if name == "main":
main()