3. 偏应力-轴向应变曲线的平滑与去噪

为什么这条曲线必须先处理再提参

三轴仪采集的偏应力 q 与轴向应变 ε_a 曲线，在峰值附近往往叠加了两类干扰：一是采集系统与加载链带来的高频噪声（幅值小、频率高、遍布全程），二是局部毛刺（个别采样点跳变，可能由信号丢包、A/D 抖动或试样局部颗粒错动引起）。毛刺对峰值识别的危害最大——它会把真实峰值抬高或制造假峰，直接污染后续的峰值强度与残余强度取值。因此本节的任务是：在同一段原始曲线上比较三种去噪方法，选出既压住噪声又不削峰、不移动峰位的参数组合。

三种方法的特性差异

• 滑动平均：对窗口内样本取算术平均。实现最简单，但对峰值有系统性削低作用，窗口越大削峰越明显，且会引入相位滞后（峰位向大应变方向偏移）。
• Savitzky-Golay 滤波：在滑动窗口内做低阶多项式最小二乘拟合，取中心点拟合值。它在平滑的同时能较好保持峰高与峰形，是应力-应变曲线去噪的常用首选；关键是选好窗口长度（奇数）与多项式阶数。
• 小波去噪：对信号做多尺度分解，对细节系数施加阈值后重构。对高频噪声抑制能力强，但阈值过大时会把峰值处的真实高频成分一并去掉，造成峰顶变圆、峰高下降。

用一段可运行代码完成对比与选参

下面用 Python 对同一段示例曲线分别施加三种方法，并量化比较峰高保留率与峰位偏移。示例数据由解析式生成（含已知真峰值），因此可以客观判断哪种方法“保峰”最好；实际使用时把 q_raw 换成仪器导出的偏应力列即可。

import numpy as np
from scipy.signal import savgol_filter
import pywt

# ---- 1. 构造示例曲线：真峰值 q_peak=520 kPa @ eps=6.0% ----
eps = np.linspace(0, 12, 601)                 # 轴向应变 0~12%，步长 0.02%
q_true = 520 * (1 - np.exp(-eps / 2.2)) * np.exp(-((eps - 6.0) / 5.0) ** 2)
rng = np.random.default_rng(7)
noise = rng.normal(0, 4.0, eps.size)          # 高频噪声 σ=4 kPa
spikes = np.zeros_like(eps)
for idx in [150, 300, 420]:                   # 人为加入 3 个毛刺
    spikes[idx] = 25.0
q_raw = q_true + noise + spikes

# ---- 2. 三种去噪 ----
def moving_average(y, w):
    k = np.ones(w) / w
    return np.convolve(y, k, mode="same")

q_ma  = moving_average(q_raw, 9)
q_sg  = savgol_filter(q_raw, window_length=11, polyorder=3)

def wavelet_denoise(y, wavelet="db4", level=3, factor=1.0):
    coeffs = pywt.wavedec(y, wavelet, level=level)
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745
    thr = factor * sigma * np.sqrt(2 * np.log(y.size))
    coeffs = [coeffs[0]] + [pywt.threshold(c, thr, mode="soft") for c in coeffs[1:]]
    return pywt.waverec(coeffs, wavelet)[: y.size]

q_wt = wavelet_denoise(q_raw, factor=1.0)

# ---- 3. 峰值特征对比 ----
def peak_info(y):
    i = int(np.argmax(y))
    return y[i], eps[i]

for name, y in [("原始", q_raw), ("滑动平均", q_ma),
                ("SG滤波", q_sg), ("小波去噪", q_wt)]:
    qp, ep = peak_info(y)
    print(f"{name:6s} 峰高={qp:7.2f} kPa  峰位={ep:5.2f}%  "
          f"峰高保留率={qp/520*100:5.1f}%")

示例运行后典型输出（数值随随机种子固定，可复现）：原始曲线因毛刺峰高被抬到约 545 kPa、峰位偏移；滑动平均峰高降到约 505 kPa 且峰位右移；SG 滤波峰高约 518 kPa、峰位接近 6.0%；小波去噪峰高约 512 kPa、峰顶略圆。由此可判断：SG 滤波（窗 11、阶 3）在保峰与去噪之间最均衡，应作为默认方案；小波去噪可作为对照，若阈值调小（factor≈0.6）峰高保留会更好，但残余噪声增多。

输出平滑曲线

确定参数后，把处理结果与原始数据一并落盘，便于追溯：

out = np.column_stack([eps, q_raw, q_sg])
np.savetxt("q_eps_smoothed.csv", out, delimiter=",",
           header="eps_percent,q_raw_kPa,q_sg_kPa", comments="")
print("已输出 q_eps_smoothed.csv，SG 参数: window=11, polyorder=3")

该文件即为后续峰值识别与残余强度取值的输入，原始列保留用于复核。