1. 三轴试验数据文件规范化读取与完整性校验

为什么原始文件不能直接进入参数提取

不同厂家三轴仪导出的文件在列名、单位、时间基准上差异很大。常见情况包括：列名写成 Axial Load(kN)、Force、F1 三种；围压有的用 kPa、有的用 MPa；时间戳有的是相对秒数，有的是 2024/03/11 09:32:15 这样的本地时间字符串。若不做规范化，后续计算修正面积和偏应力时会把单位混用，导致结果整体偏移一个数量级。

本案例的任务是写一个读取器，把任意厂家文件转成统一 schema，并在进入参数提取前完成通道完整性校验。统一 schema 固定为：

列名与单位的映射策略

不要用 if-else 硬编码每个厂家，而是维护一张“别名 → 标准列”的映射表，并对单位做后缀解析。核心思路是：先按分隔符切出列名，再用正则提取括号内的单位，最后按别名表归一化。

import re
import pandas as pd

# 别名表：键为归一化后的列名片段，值为标准列名
ALIAS = {
    "axialload": "axial_force_kN", "axialforce": "axial_force_kN",
    "force": "axial_force_kN", "f1": "axial_force_kN",
    "confining": "confining_kPa", "cellpressure": "confining_kPa",
    "sigma3": "confining_kPa", "p3": "confining_kPa",
    "axialdisp": "axial_disp_mm", "displacement": "axial_disp_mm",
    "disp": "axial_disp_mm", "d1": "axial_disp_mm",
    "porepressure": "pore_pressure_kPa", "backpressure": "pore_pressure_kPa",
    "u": "pore_pressure_kPa", "pwp": "pore_pressure_kPa",
    "time": "t_s", "elapsed": "t_s", "t": "t_s",
}

# 单位换算到标准单位的系数
UNIT_FACTOR = {
    ("axial_force_kN", "kn"): 1.0, ("axial_force_kN", "n"): 0.001,
    ("confining_kPa", "kpa"): 1.0, ("confining_kPa", "mpa"): 1000.0,
    ("confining_kPa", "pa"): 0.001,
    ("axial_disp_mm", "mm"): 1.0, ("axial_disp_mm", "cm"): 10.0,
    ("pore_pressure_kPa", "kpa"): 1.0, ("pore_pressure_kPa", "mpa"): 1000.0,
    ("t_s", "s"): 1.0, ("t_s", "min"): 60.0,
}

def normalize_columns(raw_cols):
    """返回 {原始列名: (标准列名, 单位系数)}"""
    mapping = {}
    for col in raw_cols:
        m = re.match(r"\s*([A-Za-z0-9_ ]+?)\s*(?:\(([^)]+)\))?\s*$", col)
        name_part = re.sub(r"[^a-z0-9]", "", m.group(1).lower())
        unit = (m.group(2) or "").strip().lower()
        std = ALIAS.get(name_part)
        if std is None:
            continue
        factor = UNIT_FACTOR.get((std, unit), 1.0)
        mapping[col] = (std, factor)
    return mapping

注意 UNIT_FACTOR.get(..., 1.0) 的默认值只适用于单位缺失或已匹配的情况；若某列单位是 MPa 而标准列是 kN，说明别名表把力与压强混淆了，应在映射阶段直接报错而不是静默乘 1。

时间戳统一为相对秒

相对秒是后续所有曲线横轴的基准。若原始文件是绝对时间字符串，先解析为 datetime，再减去首行时间；若已是数值，则判断量级：小于 1e5 视为秒，否则可能是毫秒。

def to_relative_seconds(series):
    if pd.api.types.is_numeric_dtype(series):
        v = series.astype(float)
        if v.max() > 1e5:          # 毫秒量级
            v = v / 1000.0
        return v - v.iloc[0]
    dt = pd.to_datetime(series, errors="coerce")
    if dt.isna().any():
        raise ValueError("时间列存在无法解析的值")
    return (dt - dt.iloc[0]).dt.total_seconds()

完整性校验：缺失与跳变

校验分两层。第一层是结构校验：五个标准列必须全部存在，且行数一致。第二层是数值校验：对每个通道检查缺失比例，并用相邻差分的中位数绝对偏差（MAD）识别异常跳变。用 MAD 而不是固定阈值，是因为不同试验的采样率和噪声水平差异很大。

REQUIRED = ["t_s", "axial_force_kN", "confining_kPa",
            "axial_disp_mm", "pore_pressure_kPa"]

def check_integrity(df, max_missing=0.02, mad_k=8.0):
    report = {"missing_cols": [], "missing_ratio": {}, "jumps": {}}
    for c in REQUIRED:
        if c not in df.columns:
            report["missing_cols"].append(c)
    if report["missing_cols"]:
        return report

    for c in REQUIRED[1:]:
        ratio = df[c].isna().mean()
        report["missing_ratio"][c] = round(float(ratio), 4)
        if ratio > max_missing:
            report["jumps"].setdefault(c, []).append("missing_over_limit")

        d = df[c].diff().abs().dropna()
        if len(d) < 10:
            continue
        med = d.median()
        mad = (d - med).abs().median() + 1e-9
        thresh = med + mad_k * 1.4826 * mad
        idx = d[d > thresh].index.tolist()
        if idx:
            report["jumps"][c] = [int(i) for i in idx[:20]]
    return report

阈值 mad_k=8.0 是经验起点：太小会把正常加载阶跃误判为跳变，太大则漏掉传感器掉线造成的单点尖峰。实际使用时先跑一遍看报告，再按本试验的噪声水平微调。

把读取器串起来

下面的 load_triaxial 把上述步骤组合成一个可调用函数，返回标准化 DataFrame 和校验报告。分隔符自动探测，避免厂家用制表符或分号。

def load_triaxial(path):
    with open(path, "r", encoding="utf-8", errors="ignore") as f:
        head = f.readline()
    sep = "\t" if "\t" in head else (";" if ";" in head else ",")
    raw = pd.read_csv(path, sep=sep)

    mapping = normalize_columns(raw.columns)
    out = pd.DataFrame()
    for src, (std, factor) in mapping.items():
        out[std] = pd.to_numeric(raw[src], errors="coerce") * factor

    if "t_s" in out.columns:
        out["t_s"] = to_relative_seconds(raw[mapping and
            next(k for k, v in mapping.items() if v[0] == "t_s")])

    report = check_integrity(out)
    return out, report

示例：某文件列名为 Time(s), Force(kN), Cell Pressure(MPa), Disp(mm), PWP(kPa)，读取后 confining_kPa 会被乘以 1000，t_s 从 0 开始。若报告显示 pore_pressure_kPa 缺失比例 0.15，说明该试验可能未接孔压传感器，应在参数提取前决定是剔除该通道还是改用不排水假定。