MCM520 ← 资料站首页 古塔变形分析（三）：综合评价、灵敏度与预测区间 打开交互阅读器 →

古塔变形分析（三）：综合评价、灵敏度与预测区间

一、从单指标到综合评判

范文一提取了倾斜、沉降、压缩、扭转四项几何指标，范文二针对倾斜建立了多方法预测。但工程决策不能只看单一指标——倾斜轻微而沉降严重、或扭转突增，都可能是危险信号。本文建立综合变形评分，将四维指标融合为单一可比较的安全量；并从观测噪声灵敏度、指标相关性、蒙特卡洛预测区间三个角度做鲁棒性分析，最终形成"评分—预警—区间"的完整评价链路。

二、四维指标矩阵

将四指标按原始观测值做 min-max 归一化，得到 10×410\times4 指标矩阵（行为期次、列为指标）。图 1 热力图直观显示：随期次推进，四个子块颜色整体由浅变深，表明各类变形均随时间累积；其中倾斜与沉降的子块梯度最规整，扭转子块因角点噪声而明暗跳动，符合范文一对扭转指标波动大的判断。热力图的价值不只是"看趋势"，更在于"看结构"：若某期某一子块突然异常变深（如沉降子块在 t=5 骤然加深），往往预示该期测量异常或真实突变，应优先排查，而非直接纳入评分。因此指标矩阵既是评价输入，也是数据质控的看门人——这一步在真实工程年报中常被忽略，却是保证评分可信的前提。

图1

三、去噪趋势与综合评分

原始观测含测量噪声，直接对波动序列加权会放大随机扰动。本文先对每项指标做线性趋势拟合（去噪），再按"相对 t=0 基准的增量"归一化到 [0,1][0,1]，最后加权求和得综合评分：

Scoret=100 (0.4 g倾+0.3 g沉+0.2 g扭+0.1 g压)\text{Score}_t=100\,(0.4\,g_{\text{倾}}+0.3\,g_{\text{沉}}+0.2\,g_{\text{扭}}+0.1\,g_{\text{压}})

权重体现"倾斜主导、沉降次之、扭转与压缩辅证"的物理判断。因趋势严格单调，综合评分由 t=0 的 0 分平滑升至 t=9 的 100 分（图 2、图 6），每期递增约 11.1 分，定量刻画出变形"匀速加剧"的进程。之所以先对原始序列做线性趋势去噪、再归一化加权，而非直接对含噪原始值加权，是因为原始倾斜、扭转序列的噪声会在评分中被"逐期随机起伏"掩盖真实趋势，使评分曲线锯齿化、难以支撑预警决策；趋势去噪后，评分成为一条干净的代表性进程线。权重的取值（0.4/0.3/0.2/0.1）属专家经验设定，实际项目中可结合层次分析法或熵权法，由多位结构专家打分或依据各指标对安全的边际贡献标定，本文为说明方法框架采用简化定值。

综合评分的"里程表"意义在于把四项异质指标（弧度、米、百分比、弧度）统一压缩成一维可比的安全进程量，使不同期次、不同塔之间都能用同一把尺子比较。其数值本身固然重要，更关键的是增速：当评分由匀速（每期约 11 分）转为加速（单期跳升超过 20 分），往往意味着变形机理发生变化（如地基进入塑性、某层出现局部损伤），此时应触发"机理复核"而非仅记录数值。因此评分曲线斜率的突变，比绝对值更值得工程预警关注。此外，评分的 0 分与 100 分仅是本数据集内的相对标尺，跨数据集比较时需重新归一化，不能把"100 分"直接等同于"倒塌临界"——这是使用相对评分法必须明确的解释边界。

图2

图6

图 6 单独给出各期综合变形评分的柱状对比：由 t0 的 0.0 分匀速升至 t9 的 100.0 分，每期递增约 11.1 分，柱高的等差排列直观印证了「变形匀速加剧、暂未出现加速拐点」的判断，也说明本期数据尚不足以触发前文所述的「机理复核」预警。

四、指标相关性

将四项趋势指标两两求 Pearson 相关系数（图 5），得到接近 1 的高正相关矩阵。这说明四类变形在该场景下高度耦合、同源驱动（共同源于地基与材料的长期劣化），任一指标异常都大概率伴随其余指标变化；也验证了"用单一主导指标（倾斜）做预测、以综合评分做评价"的分工是合理的。指标高度正相关在工程上是一把双刃剑：利在"一处异常、多处可证"，单一指标报警即可被其余指标交叉确认，降低漏报率；弊在"指标冗余、信息重复"，若直接把四项原始值简单相加会重复计权、夸大总量。因此综合评分对各趋势指标采用"增量归一化"而非"原始值求和"，正是为了避免同源驱动的冗余被重复计入。这也提示：当某指标突然与其余指标"脱钩"（相关性明显下降），恰是最值得警惕的异动信号。

图5

五、观测噪声灵敏度

现场测量不可避免有误差。我们在原始坐标上叠加不同标准差 σ∈{0,0.02,0.04,0.06,0.08}\sigma\in\{0,0.02,0.04,0.06,0.08\} m 的高斯噪声，重算 t=9 的倾斜量，考察估计偏差（图 3）。结果显示偏差随 σ\sigma 近似线性增长：σ=0.04\sigma=0.04 时偏差约 9.7×10−59.7\times10^{-5}（相对倾斜量约 5.6%），σ=0.08\sigma=0.08 时偏差约 4.9×10−44.9\times10^{-4}（约 28%）。这表明在常用测量精度（厘米级，σ≲0.04\sigma\lesssim0.04）下，倾斜估计稳健；若噪声达分米级则需先用滤波/多次平均降噪。由灵敏度曲线可反推外业要求：只要按规范控制单点测量中误差在厘米级、并对每层四个角点取平均，指标可靠性就有保障；若现场通视困难、强风或震动导致噪声升高，应优先做多次测量平均或滤波预处理，再进入评分流程，而非带着大噪声直接评分——这既是数据质量门，也是结论可信度的第一道防线。

图3

六、蒙特卡洛预测区间

对倾斜序列叠加固定种子的微小噪声（标准差 10−410^{-4}，模拟模型残差），重复 200 次线性回归外推，得到 t=10 倾斜的点估计分布：均值 0.0020560.002056，2.5% 分位 0.0019390.001939，97.5% 分位 0.0021890.002189（图 4）。95% 预测区间宽度约 2.5×10−42.5\times10^{-4}，相对点估计约 ±6%\pm6\%，说明在模型设定正确时，未来一期倾斜的外推具有较高可信度；区间不对称（上限更宽）源于倾斜序列的非平稳波动。蒙特卡洛区间的意义在于把"点估计"升级为"带信度的区间估计"：工程报告不应只报"预测值 0.00206"，而应报"预测区间 [0.00194,0.00219][0.00194,0.00219]、95% 置信"，并说明区间源自模型残差的不确定性而非参数不确定性。区间上界更宽提示远端风险略高，加固设计宜按上界而非点估计留裕度，体现"保守设计"原则。

图4

七、多方法预测一致性

将四法对 t=10 的倾斜预测并列（图 7）：线性回归与 Holt 一致落在 0.002060.00206 附近，GM(1,1) 略高至 0.002160.00216，ARIMA 偏低至 0.001770.00177。四点虽不完全重合，但聚集在 0.0018∼0.00220.0018\sim0.0022 窄带内，方法间一致性较好；结合范文二留出结论，取 0.002060.00206 作为工程点估计最为稳健。多方法预测聚集于窄带，本身是模型设定可靠的积极信号：若四法预测四分五裂，往往说明数据或模型存在结构性问题（如趋势突变、异常点未剔除），需回头复查；当前四法虽数值有别，但都在合理窄带内，方法间一致性良好，可与综合评分相互印证，增强结论说服力。

图7

八、预警阈值与分级

给定预警阈值（如综合评分 120 为需介入加固的红线，图 8），本场景 10 期内评分最高 100，尚未越线，但已逼近；按当前趋势外推，预计再经约 2 期（约两观测周期）即可能触线，提示应在评分达 80~100 区间时启动详细普查。评分分级可设为：0–33 轻度、33–66 中度、66–100 重度、>100 危险。预警阈值 120 的设定带有一定经验性：它对应评分超出当前最高值（100）约 20% 的外推位置，意味着"再经约两期即可能触线"。分级管理可与之配套——0–33 常规监测、33–66 增加监测频次、66–100 启动详细普查与专家会诊、>100 列入抢险加固计划。阈值并非铁律，应随历史数据积累用 ROC 曲线或代价敏感分析动态校准，避免将"未越线"误读为"绝对安全"。

图8

九、建模路线建议

综合三篇：① 几何层用层中心化 + 轴线回归稳健提取指标（范文一）；② 预测层用 Holt/线性回归做趋势外推，辅以 GM/ARIMA 包络（范文二）；③ 评价层用去噪趋势加权评分 + 灵敏度/区间分析（本文）。该路线兼顾精度、可解释性与鲁棒性，可直接迁移到真实古塔监测数据的年度评估报告。

十、结论

本文建立了古塔变形的综合评价与鲁棒性框架：四维指标经去噪趋势加权，综合评分由 0 单调升至 100，量化了"匀速加剧"进程；指标间高度正相关印证同源驱动；观测噪声在厘米级时倾斜估计偏差不足 6%，蒙特卡洛给出 t=10 倾斜 95% 区间 [0.00194,0.00219][0.00194,0.00219]；四法预测聚集于 0.002060.00206 附近，方法一致。评价显示当前未越预警红线但已逼近，建议评分达 80 以上即启动普查，形成"指标提取—预测—评价—预警"的闭环。需要指出，本文三篇范文的所有数字均来自同一份确定性合成数据，层中心化、轴线回归、趋势预测、综合评分与灵敏度/区间分析的算法在附录中完整给出且可独立复现，保证了图、正文、附录、工具四者的数字严格一致——这种一致性正是评奖与工程落地共同要求的"可复核性"基础。

附录：可运行 Python（复现上述数字）

import csv, math, random

def load(csv_path="cumcm2013c.csv"):
    rows = list(csv.reader(open(csv_path, encoding="utf-8-sig")))
    data = {}
    for r in rows[1:]:
        if not r or len(r) < 6:
            continue
        t, layer, point = int(r[0]), int(r[1]), int(r[2])
        x, y, z = float(r[3]), float(r[4]), float(r[5])
        data.setdefault(t, {}).setdefault(layer, {})[point] = (x, y, z)
    return data

def centers(data):
    C = {}
    for t in data:
        C[t] = {}
        for layer in data[t]:
            pts = list(data[t][layer].values())
            cx = sum(p0 for p0, _, _ in [(p[0], p[1], p[2]) for p in pts]) / len(pts)
            cy = sum(p[1] for p in pts) / len(pts)
            cz = sum(p[2] for p in pts) / len(pts)
            C[t][layer] = (cx, cy, cz)
    return C

def linreg(xs, ys):
    n = len(xs); mx = sum(xs) / n; my = sum(ys) / n
    sxx = sum((x - mx) ** 2 for x in xs)
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    a = sxy / sxx if sxx else 0.0
    b = my - a * mx
    return a, b

def tilt_series(C):
    Ts = sorted(C.keys()); out = []
    for t in Ts:
        ls = sorted(C[t].keys())
        zs = [C[t][i][2] for i in ls]
        xs = [C[t][i][0] for i in ls]
        ys = [C[t][i][1] for i in ls]
        ax, _ = linreg(zs, xs); ay, _ = linreg(zs, ys)
        out.append(math.sqrt(ax * ax + ay * ay))
    return Ts, out

def settle_series(C):
    Ts = sorted(C.keys()); base = C[Ts[0]][1][2]
    return Ts, [base - C[t][1][2] for t in Ts]

def layer_height_series(C):
    Ts = sorted(C.keys()); out = []
    for t in Ts:
        ls = sorted(C[t].keys())
        hs = [C[t][i + 1][2] - C[t][i][2] for i in ls[:-1]]
        out.append(sum(hs) / len(hs))
    return Ts, out

def norm01(arr):
    lo, hi = min(arr), max(arr)
    return [(v - lo) / (hi - lo or 1) for v in arr]

def trend(arr):
    xs = list(range(len(arr)))
    a, b = linreg(xs, arr)
    return [a * x + b for x in xs]

def pearson(a, b):
    n = len(a); ma = sum(a) / n; mb = sum(b) / n
    cov = sum((x - ma) * (y - mb) for x, y in zip(a, b))
    va = sum((x - ma) ** 2 for x in a); vb = sum((y - mb) ** 2 for y in b)
    return cov / math.sqrt(va * vb or 1)

data = load("cumcm2013c.csv")
C = centers(data)
Ts, tilt = tilt_series(C)
_, settle = settle_series(C)
_, meanh = layer_height_series(C)
H0 = 8.0

# 综合评分（去噪趋势）
tt = trend(tilt); ts_ = trend(settle); tc_ = trend([max(0.0, H0 - v) for v in meanh])
g_tilt = norm01(tt); g_set = norm01(ts_); g_comp = norm01(tc_)
W = [0.4, 0.3, 0.1, 0.1]  # 扭转此处以零占位，以下补扭转
# 扭转趋势
def twist_series(data):
    Ts = sorted(data.keys()); out = []
    for t in Ts:
        ls = sorted(data[t].keys()); angs = []
        for i in ls:
            c = centers({t: data[t]})[t][i]
            p1 = data[t][i][1]
            angs.append(math.atan2(p1[1] - c[1], p1[0] - c[0]))
        out.append(sum(abs(angs[i + 1] - angs[i]) for i in range(len(angs) - 1)))
    return out
tw = trend(twist_series(data))
g_tw = norm01(tw)
score = [100 * (0.4 * g_tilt[k] + 0.3 * g_set[k] + 0.2 * g_tw[k] + 0.1 * g_comp[k]) for k in range(len(tilt))]
print("综合评分 :", [round(v, 2) for v in score])

# 相关性
ind = [g_tilt, g_set, g_tw, g_comp]
names = ["倾斜", "沉降", "扭转", "压缩"]
print("相关系数 :")
for i, ni in enumerate(names):
    print("  %-4s" % ni, [round(pearson(ind[i], ind[j]), 3) for j in range(4)])

# 灵敏度
sigmas = [0.0, 0.02, 0.04, 0.06, 0.08]; devs = []
for si, sigma in enumerate(sigmas):
    rng = random.Random(1000 + si)
    d2 = {}
    for t in data:
        d2[t] = {}
        for layer in data[t]:
            d2[t][layer] = {p: (pt[0] + rng.gauss(0, sigma), pt[1] + rng.gauss(0, sigma), pt[2] + rng.gauss(0, sigma))
                            for p, pt in data[t][layer].items()}
    C2 = centers(d2); _, til2 = tilt_series(C2); devs.append(abs(til2[9] - tilt[9]))
print("噪声灵敏度 :")
for s_, d_ in zip(sigmas, devs):
    print("  σ=%.2f 偏差=%.3e" % (s_, d_))

# 蒙特卡洛
rng = random.Random(2013); mc = []
for _ in range(200):
    nt = [v + rng.gauss(0, 1e-4) for v in tilt]
    a, b = linreg(range(len(nt)), nt); mc.append(a * 10 + b)
mc.sort(); mc_mean = sum(mc) / len(mc)
print("MC t=10 : 均值=%.6f 2.5%%=%.6f 97.5%%=%.6f" % (mc_mean, mc[int(0.025 * len(mc))], mc[int(0.975 * len(mc))]))