企业官网建设流程全解析

简介：分段多项式拟合在动三轴试验数据处理中的应用是一篇PDF格式的学术论文，聚焦海洋平台地基土动三轴试验中离散、有限数据点的拟合难题，适合从事土动力学、岩土工程数据处理的研究人员与工程师阅读参考。文件包共1份PDF，大小836KB，内容涵盖最小二乘拟合原理、分段低次多项式方法、剪应变压缩处理及正规方程病态问题的改进思路。基于渤海某海洋平台地基土试验数据，作者采用该方法拟合动剪模量、阻尼比与剪应变的关系，并与Hardin-Black、黄雨、陈国兴等经验公式对比，论证了分段多项式拟合的适用性与可靠性，为获取不同剪应变下土动力特性参数提供了有价值的计算参考。已有123人浏览/学习，可结合原文公式与图表对照学习试验数据拟合全过程。

1. 分段多项式拟合在动三轴试验数据处理中的角色

分段多项式拟合处理动三轴试验数据，第一步不是写拟合函数，而是分清滞回圈里的加载段和卸载段。循环应力-应变轨迹是闭合曲线，同一应变处加载与卸载对应不同应力，不是单值函数，一条高次多项式整体拟合必然失败。

分段多项式拟合按物理转折点把曲线切开，每段用低次多项式独立拟合，段间加连续性约束，再用拟合曲线计算动弹性模量与阻尼比。用 NumPy 与 SciPy 就能实现，不依赖商用岩土软件。

下文按转折点识别、分段拟合、参数提取、自适应精修的思路展开，代码直接针对动三轴滞回圈数据，参数建议来自实际处理经验。

2. 动三轴滞回圈的转折点识别：分段位置的确定

2.1 动三轴试验数据为什么不能整体拟合

动三轴试验在恒定围压下对试样施加循环轴向应力 σd，同步记录轴向应变 εd。一个完整循环的 σd-εd 轨迹就是滞回圈：应力从谷值单调升到峰值，再单调回落到谷值。滞回圈面积代表一个循环中土体损耗的能量，是阻尼比计算的分子项，所以数据处理的第一步是把加载段和卸载段准确分开，也就是找到转折点。

物理上，加载段是土体继续压缩、骨架刚度不断调整的过程；卸载段是弹性回弹与塑性恢复并存的过程。加载和卸载在重叠的应变区间内对应不同应力，σd 与 εd 不是单值函数关系，用一条整段多项式去表达两个分支，无论次数多高都得不到同时贴合两段的解。

数值上，对整条曲线做高次多项式拟合，在数据区间端点附近会出现明显摆动。动模量恰好要用峰、谷两端的应力应变差来算，端点一翘，模量和阻尼比全部偏掉。分段多项式拟合把问题拆成两个子问题：节点位置在哪、每段用什么多项式，下文先解决节点位置。

2.2 用应力差分变号定位转折点

转折点就是应力时间序列的极值点：应力对时间的一阶差分由正变负处是峰值，由负变正处是谷值。动三轴采样率常见 100 Hz～1 kHz，直接差分会被传感器噪声打出一串假极值，因此先做移动平均，再用最小间隔过滤候选点。

import numpy as np def find_turning_points(stress, min_gap=10, smooth_window=3): """ 通过应力差分符号变化定位滞回圈转折点 stress : 一次循环的应力时间序列，一维数组 min_gap : 相邻转折点之间的最少采样点数，用于滤除噪声 smooth_window: 差分前移动平均窗口，建议 3~7，过大会抹平真实极值 返回转折点在原始数组中的索引，按时间顺序排列 """ if smooth_window > 1: kernel = np.ones(smooth_window) / smooth_window sig = np.convolve(stress, kernel, mode='same') else: sig = stress dsigma = np.diff(sig) sign = np.sign(dsigma) sign_change = np.where(np.diff(sign) != 0)[0] + 1 valid = [] for idx in sign_change: if not valid or idx - valid[-1] >= min_gap: valid.append(idx) return np.array(valid)

逻辑说明：np.convolve 用三点平均窗压掉高频抖动；np.diff 计算一阶差分，np.sign 映射成 -1/0/1；np.diff(sign) != 0 出现的位置就是符号翻转点，加 1 是因为 diff 结果比原数组短一位。smooth_window 取奇数，mode='same' 下输出与输入索引对齐，不会引入半个窗口的偏移。

min_gap 按采样率调整：一个循环 500 点时取 10~20，1000 点时取 30~50。若差分出现平台（应力连续多个采样点不变），平台两个边缘都会触发变号，min_gap 会把它们过滤成一组，只保留平台起点。波形畸变严重时，用应力极值切分会把相位差带来的回折段切进分支，常见做法是把入参换成 strain，找 dε/dt 变号点；两个结果对比差值超过 5 个采样点时，说明试样粘滞明显，以应变变号点为准。

2.3 转折点物理含义与滞回圈分割约定

两个转折点把一次循环分成两段：应力峰值之前的上升段是加载段，峰值之后的下降段是卸载段。工程上习惯用应力极值定义加载/卸载，对应关系如下表。

多级动三轴试验一个试样几十个循环，转折点还承担圈与圈的分割：stress[转折点[i] : 转折点[i+1]] 就是一段完整分支。首圈常带初始压实段，末圈可能不闭合，进入拟合前把首尾各 10 个点裁掉，避免把开口偏差带进面积积分。

3. 截断幂基实现分段多项式拟合：最小二乘与样条

3.1 截断幂基与段间连续性

分支内部如果只用一个多项式，加载段这种“低应变高刚度、高应变低刚度”的双线性形态拟合不好，常见做法是分支内部再分段，节点取刚度拐点对应的应变。固定节点位置的分段多项式拟合用截断幂基：设节点为 k₁ 到 k_m，段内次数 D，模型写成 y = β₀ + Σβ_d·x^d + ΣᵢΣ_d γᵢ,ᵈ·(x−kᵢ)+^d，其中 (u)+ = max(u, 0)。

这个基函数天然满足段间 C0 连续：每个截断项在节点处取值 0，左右两段曲线在节点上必然相交，不需要额外写边界条件。D=2 时节点两侧斜率有折点，正适合滞回圈这种“加载段与卸载段斜率不同”的形态；D=3 再追加 (x−k)_+³ 项，一阶导也连续，适合带回弹弧线的卸载分支。

所有系数进入同一个线性最小二乘问题，不需要分段循环求系数，也不会出现分段接缝处的振荡拼装。参数个数 p = 1 + D + D·m，动三轴一个分支通常几百到上千个采样点，自由度充裕。

3.2 固定节点的最小二乘拟合代码

def fit_piecewise_poly(x, y, knots, degree=2): """ 截断幂基分段多项式拟合，自动满足段间 C0 连续 x, y : 一条滞回圈分支的应变与应力 knots : 内部节点位置的 x 坐标，需递增且在数据区间内部 degree : 每段多项式次数，2 或 3 返回系数 coef 与可调用函数 predict """ cols = [np.ones_like(x)] for d in range(1, degree + 1): cols.append(np.power(x, d)) for k in knots: for d in range(1, degree + 1): cols.append(np.maximum(x - k, 0.0) ** d) A = np.column_stack(cols) coef, *_ = np.linalg.lstsq(A, y, rcond=None) def predict(xx): xx = np.asarray(xx, dtype=float) val = np.full_like(xx, coef[0], dtype=float) for d in range(1, degree + 1): val += coef[d] * np.power(xx, d) ptr = degree + 1 for k in knots: for d in range(1, degree + 1): val += coef[ptr] * np.maximum(xx - k, 0.0) ** d ptr += 1 return val return coef, predict

设计矩阵 A 的第一列是常数项，随后是 x 到 x^D；对每个节点 k 追加 (x−k)+ 到 (x−k)+^D 共 D 列。np.linalg.lstsq 解最小二乘，predict 闭包按相同基重建曲线，供第 4 章取端点值、算面积。knots 传应变坐标，第 2.2 节返回的是索引，用 strain[idx] 转换。

knots 不能落在数据区间端点：节点与端点重合时截断项在区间内恒为 0，设计矩阵出现列相关。赋值前用 np.clip(knots, x.min()+Δ, x.max()−Δ) 约束，Δ 取应变幅值的 1% 即可。knots 有多个时要递增，否则基函数顺序错乱。

3.3 用 UnivariateSpline 做分支平滑拟合

节点位置拿不准或不关心物理含义时，常见做法是对每个分支单独做三次样条平滑，节点由算法自动选。

from scipy.interpolate import UnivariateSpline def fit_branch_spline(strain, stress, s=None): """s 为平滑因子，量纲是残差平方和；s=0 时穿过全部数据点，不推荐""" n = len(strain) if s is None: s = n * np.var(stress) * 0.05 spl = UnivariateSpline(strain, stress, k=3, s=s) return spl

s 设太小，样条追随传感器噪声，转折点附近出现小抖动；设太大，滞回圈被抹平成一条线，阻尼比明显偏小。经验区间是 n×var(stress) 的 2%~10%，选完后用面积偏差校验（见第 5 章）。样条适合只求曲线平滑的场景；报告里需要写明“分段位置”时，用 3.2 节的固定节点方式，物理含义更明确。

3.4 段内次数与连续性选型表

次数不是越高越好，D=4 时截断项对节点邻域误差敏感，节点两侧出现波浪形过拟合。滞回圈的物理特征是斜率在转折点突变，D=2 或 3 已经把特征表达完整，再高只会放大噪声。

4. 从分段拟合结果提取动弹性模量与阻尼比

4.1 割线动模量的计算口径

动弹性模量反映循环荷载下的刚度，最常用的是割线模量：峰值应力减谷值应力，除以峰值应变减谷值应变。峰值应力取拟合曲线在转折点处的函数值而不是原始测点值，单点噪声被过滤，这是分段拟合的第一个落点。

def dynamic_modulus(eps_peak, sig_peak, eps_valley, sig_valley, strain_in_percent=True): """割线动模量 Ed = (σ峰值 − σ谷值) / (ε峰值 − ε谷值) 输入建议取 3.2 节 predict 或 np.polyval 在转折点处的拟合值 strain_in_percent=True 时先换算回小数""" if strain_in_percent: eps_peak /= 100.0 eps_valley /= 100.0 return (sig_peak - sig_valley) / (eps_peak - eps_valley)

单位统一为 kPa：应力传感器给 MPa 就先乘 1000，应变给 % 就除 100，模量换算错误是结果表里最常见的错误。切线模量由拟合多项式解析求导得到，取加载段起始区间的 dσ/dε，用于小应变幅值（10⁻⁴ 以下）的初始刚度估计；常规应变幅值下报告用割线模量。

4.2 滞回圈面积与阻尼比公式

阻尼比 λ = A_loop / (4π·W)，A_loop 是滞回圈面积，W 是半循环内储存的弹性能，对线性弹性段 W = 0.5·σa·εa，σa、εa 是应力、应变幅值。面积从拟合分支做差积分：A_loop = ∫(加载分支 − 卸载分支) dε，用高密度采样梯形积分实现。

def hysteresis_loop_area(eps_min, eps_max, load_predict, unload_predict): """滞回圈面积 = |∫(加载分支 − 卸载分支) dε|，应变方向统一递增""" x = np.linspace(eps_min, eps_max, 5000) return np.abs(np.trapz(load_predict(x) - unload_predict(x), x)) def damping_ratio(loop_area, sig_amp, eps_amp): """阻尼比 = 滞回圈面积 / (4π × 弹性能)，σa、εa 为幅值""" W = 0.5 * sig_amp * eps_amp if W <= 0: return np.nan return loop_area / (4.0 * np.pi * W)

注意三个坑。第一，积分方向用“加载减卸载”，写反了面积变负，取 abs 兜底。第二，弹性能用的是幅值，如果手头只有峰谷差 Δσ、Δε，等效式是 λ = 2·A_loop/(π·Δσ·Δε)，两个式子结果必须一致，正好用来核对代码。第三，原始测点在峰值附近有毛刺时，梯形积分会放大误差，拟合分支在节点处受最小二乘约束，面积对单点噪声不敏感。新版 NumPy 2.x 把 np.trapz 改名为 np.trapezoid，报错时替换即可。

4.3 多级动三轴数据的批量处理

多级试验文件结构通常是级别编号、时间、应力、应变。批量处理按级别分组，每级内部找转折点，每两个相邻转折点切出一段分支，用 np.polyfit 做二次多项式拟合：

import pandas as pd def batch_process(csv_path, min_gap=10, degree=2): df = pd.read_csv(csv_path) rows = [] for level, g in df.groupby('level'): strain = g['strain'].values stress = g['stress'].values idx = find_turning_points(stress, min_gap=min_gap) if len(idx) < 4: # 至少要有峰值、谷值、下一个峰值 continue for i in range(1, len(idx) - 1): x = strain[idx[i]:idx[i + 1]] y = stress[idx[i]:idx[i + 1]] coef = np.polyfit(x, y, degree) # 每支独立二次拟合 y_fit = np.polyval(coef, x) r2 = 1 - np.sum((y - y_fit) ** 2) / np.sum((y - y.mean()) ** 2) # 按应力升降方向判断分支起止，避免谷峰颠倒 if stress[idx[i]] < stress[idx[i + 1]]: # 上升段：谷到峰 eps_valley, eps_peak = strain[idx[i]], strain[idx[i + 1]] else: # 下降段：峰到谷 eps_valley, eps_peak = strain[idx[i + 1]], strain[idx[i]] rows.append({ 'level': level, 'branch_index': i, 'eps_peak': eps_peak, 'eps_valley': eps_valley, 'sig_peak': np.polyval(coef, eps_peak), 'sig_valley': np.polyval(coef, eps_valley), 'r_squared': r2, }) return pd.DataFrame(rows)

分支形态是简单的单调段时，单支二次多项式足够，不必上截断幂基；只有分支内部出现明显拐点才需要 3.2 节的内部节点。加载段和卸载段不要求完全对称，逐支拟合保留各自形态。加载分支本身从谷值到峰值，单条记录就够算 Ed；阻尼比需要把相邻的加载、卸载两条分支成对带入 4.2 节的两个函数。结果字段约定见下表，落库时统一 long format，方便画模量衰减曲线。

5. 分段多项式拟合的段数与节点自适应精修

5.1 用 AIC 选段内多项式次数

段内次数默认 2，数据形态偏离时用 AIC 在 1~4 之间选。AIC 在残差与参数个数之间折中：

def select_degree_aic(x, y, knots, max_degree=4): """比较不同段内次数的 AIC，返回 (aic, 最优次数, 参数个数, rss)""" n = len(x) best = None for deg in range(1, max_degree + 1): coef, predict = fit_piecewise_poly(x, y, list(knots), degree=deg) rss = float(np.sum((predict(x) - y) ** 2)) p = len(coef) aic = n * np.log(rss / n) + 2 * p if best is None or aic < best[0]: best = (aic, deg, p, rss) return best

AIC 越小越好。节点多、次数高时 p 迅速增大，惩罚项压住高次模型；噪声大的数据常选回 degree=1，说明节点划分已经逼近噪声下限。分支内部要判断“是否再多切一段”时，把候选节点按应变分位数（25%、50%、75%）分别带入 select_degree_aic，AIC 最小的那一组同时给出节点取舍和段内次数。

5.2 转折点位置的亚采样精修

2.2 节差分法只定位到采样点索引，真实应力峰值往往落在两个采样点之间。用局部抛物线拟合把转折点精修到小数索引：

def refine_turning_index(stress, idx, half_window=5): """在初步转折索引附近拟合局部抛物线，返回亚采样精度的峰值位置""" lo = max(0, idx - half_window) hi = min(len(stress), idx + half_window + 1) t = np.arange(lo, hi) - idx coef = np.polyfit(t, stress[lo:hi], 2) apex = -coef[1] / (2.0 * coef[0]) if abs(apex) > half_window: # 顶点跑出窗口说明不是单峰 return float(idx) return idx + apex

抛物线顶点就是局部极值的精确位置。abs(apex) > half_window 说明窗口内不是平滑单峰（例如应力平台或毛刺），保持原索引。精修之后是小数索引，取应变节点时用 np.interp(refined, np.arange(n), strain)，再把该值传给 3.2 节 knots。这一处修正对阻尼比影响明显：峰值点偏 2~3 个采样点，面积积分偏差会在 1% 量级。

5.3 拟合质量的三个校验判据

拟合完先看 R² 之外的三件事：

原始梯形积分面积用 np.trapz(stress, strain) 预计算，与 4.2 节拟合面积对比。偏差持续大于 1% 时优先怀疑节点位置而不是段内次数：把 5.2 节精修后的节点替换回去重跑。残差序列只在节点处出现 V 形尖峰，说明该处斜率突变是模型强加的，考虑在该节点追加 (x−k)³ 项换成平滑过渡。

6. 动三轴数据批量落地值得留意的 3 个细节

6.1 按循环编号切圈，不跨圈拟合

多级试验一个文件几十个循环，整体找转折点会把不同循环间的残余位移算进同一分支。处理前用 scipy.signal.find_peaks 在应力时间序列上按周期找峰，峰与峰之间是一圈，每圈单独执行 2.2 与 3.2 节流程。首圈初始压实段、末圈不闭合段裁掉再拟合。

6.2 用 sidecar 文件记录处理参数

min_gap、degree、s、节点位置这些参数直接影响模量与阻尼比。处理后在同一目录输出 JSON 后缀文件，记录数据文件名、转折点数、段内次数、AIC 值、面积偏差。实验室间比对时，能复现处理过程的结果才有讨论价值。

6.3 面积积分优先用拟合曲线

滞回圈开口会让 np.trapz(stress, strain) 产生虚假面积。拟合积分以转折点为边界，天然忽略开口偏移，但对转折点检测失败敏感。批量脚本里把面积偏差阈值设为 1%，超限数据不进结果表，单独导出供人工复核。

最后留一个可落地的细节：批量输出前把每圈的拟合参数、残差均值、面积偏差写进同一份 CSV，调试时按偏差降序排列，优先看偏差最大的几条曲线，通常能找到传感器滑零或增益漂移。这个顺序比按圈号顺序检查省时间。

本文还有配套的精品资源，点击获取