古塔变形分析（二）：倾斜趋势的时间序列预测

一、从指标到预测

范文一从几何上提取了倾斜、沉降、压缩、扭转四项变形指标。对古塔保护最有决策价值的是趋势外推：已知前 10 期（t=0…9）的变形指标，如何预测未来若干期（t=10,11,12…）的变形量，从而判断何时达到预警阈值、安排加固。

本文以轴线倾斜量 $\tan\theta$ 为核心预测对象（它综合反映塔体侧倾，且趋势最稳健），系统比较四种时间序列预测方法：线性回归、GM(1,1) 灰色预测、Holt 指数平滑、ARIMA(1,1,0)，并用留出验证（留最后两期做检验）评估精度，最终给出未来三期的预测区间。

二、预测对象与数据形态

在正式建模前，有必要明确"预测什么"与"为何预测"。古塔监测的终极目标是判断"何时需要干预"，而这依赖于对变形速率的量化外推。倾斜量之所以成为首选预测目标，除前文所述敏感性、稳健性、耦合性三重理由外，还因为它具有清晰的物理单位（无量纲的 tanθ 或可直观换算的顶层偏移米数），便于与非专业的文保、监理人员沟通；相比之下，扭转 rad 或压缩率 % 的预警阈值更难被现场人员直观理解。因此，以倾斜量为主预测量、辅以其余指标校核，既是数学上的优解，也是工程沟通上的务实选择。下文四种方法均围绕该序列展开。

倾斜量序列（由范文一层中心回归得到）为：
$\tan\theta = [0.001038, 0.001164, 0.000836, 0.001238, 0.001531, 0.001345, 0.001508, 0.002145, 0.001868, 0.001745]$ 。

该序列含测量噪声、整体上升，属典型的"带噪趋势型"时间序列。图 1 绘出原始序列与线性回归拟合线（虚线外推至 t=12），可见线性模型已能捕捉主体趋势，斜率 $0.0001119$ /期。选择倾斜量而非沉降或扭转作为预测主角，有三重理由：其一，倾斜是评价古塔稳定性的最敏感指标，顶层偏移对其一阶敏感；其二，倾斜序列趋势最稳健、噪声相对最小，模型外推可信度最高；其三，倾斜与沉降、扭转高度耦合，预测倾斜即可间接反映整体变形走向。这也符合工程监测"抓主导、带全局"的务实原则——不必对四个指标各建一套预测，而应以主导指标牵引、其余指标校核。

三、方法一：线性回归

将期次 $t$ 作为自变量，对 $\tan\theta_t=a\,t+b$ 做最小二乘，得 $a=0.00011187,\ b=0.00093836$ 。外推：
$\hat{\tan\theta}_{10}=0.002057,\ \hat{\tan\theta}_{11}=0.002169,\ \hat{\tan\theta}_{12}=0.002281$ 。

线性模型假设变形匀速发展，简洁透明、可解释性强，是工程预测的常用基线。其最大优势在于参数极少（仅斜率与截距）、运算完全透明，能直接给出工程最关心的"每期增量"——本文为 $1.12\times10^{-4}$ /期；劣势则是对非线性加速不敏感，若后期变形因某种诱因突然加快，线性会系统性低估。因此线性模型最适合"中短期、趋势平稳"的外推，并作为基线与其他非线性方法对照，而非独立决策依据。

四、方法二：GM(1,1) 灰色预测

灰色预测适用于"少数据、贫信息"场景。对原始序列做一次累加生成（AGO），在生成序列上建立一阶微分方程 $\mathrm{d}x^{(1)}/\mathrm{d}t+a x^{(1)}=u$ ，用最小二乘估计 $(a,u)$ ，再累减还原。图 3 显示 GM(1,1) 拟合与原序列吻合良好，其外推为：
$\hat{\tan\theta}_{10}=0.002158,\ \hat{\tan\theta}_{11}=0.002334,\ \hat{\tan\theta}_{12}=0.002524$ 。

GM(1,1) 因指数结构，外推增速略快于线性，反映其捕捉到的"加速"成分。灰色预测的理论动机在于：它假定原始序列经一次累加后近似指数规律，从而可用一阶微分方程描述。对单调增长且样本稀少的数据，它能比线性更好地拟合"轻微加速"；但其弱点恰在多步外推——指数函数在远端急剧放大，一旦真实趋势并非严格指数（如受修复工程干预而减速），预测将严重偏离。本文留出验证中 GM 误差最大（MAPE 28.4%），正是这一弱点的体现：它把噪声扰动误读为"加速信号"并过度外推。故 GM(1,1) 宜用于"样本内拟合优度"展示，外推决策须谨慎。

五、方法三：Holt 指数平滑

Holt 双参数法同时建模水平与趋势：对每个时刻更新 level 与 trend，预测值为 level + trend·h。取 $\alpha=0.3,\ \beta=0.1$ ，拟合曲线如图 4，外推为：
$\hat{\tan\theta}_{10}=0.002059,\ \hat{\tan\theta}_{11}=0.002171,\ \hat{\tan\theta}_{12}=0.002282$ 。

Holt 结果几乎与线性回归一致，说明趋势项主导，平滑仅微调近期权重。Holt 可视为"带趋势项的指数平滑"：水平分量跟踪当前值，趋势分量跟踪变化率，二者通过 $\alpha,\beta$ 平衡历史与近期信息。它比纯线性多了一个"趋势自适应"自由度，又比 GM 少了强指数假设，因而在"有趋势、有噪声、样本少"场景下往往最稳健——本文留出验证中它 MAPE 最低（10.8%）即印证。 $\alpha=0.3,\beta=0.1$ 为经验默认，实际可按留一出交叉验证进一步微调以适配具体数据。

六、方法四：ARIMA(1,1,0)

先对序列作一阶差分消除趋势，再对差分序列建 AR(1) 模型 $d_t=\phi d_{t-1}+e$ 。差分后序列近似平稳，估计得 $\phi\approx0$ （噪声主导），故 ARIMA 预测趋于"回归到末值附近"：
$\hat{\tan\theta}_{10}=0.001769,\ \hat{\tan\theta}_{11}=0.001764,\ \hat{\tan\theta}_{12}=0.001765$ 。

图 5 显示 ARIMA 拟合在前段贴近、末段保守，外推值明显低于其余三法，提示其对"趋势持续性"假设较弱。ARIMA(1,1,0) 先差分消除趋势、再对平稳差分序列建自回归，是经典 Box-Jenkins 范式。它对平稳或周期序列极为有效，但面对"强趋势、短样本"的古塔数据，一阶差分几乎抹去趋势信息，残差近似白噪声，导致外推退化为"延续末值"，故预测偏保守（约 0.00177）。若要发挥 ARIMA 优势，应升阶（如 ARIMA(1,2,0) 二次差分或加入季节项），或先对原序列去趋势再建模；在样本仅十余期的古塔场景，其性价比通常不及 Holt 与线性回归。

七、留出验证与精度对比

为客观比较，用前 8 期（t=0…7）训练，预测 t=8、t=9 两个实际值，计算 RMSE 与 MAPE：

方法	RMSE	MAPE
线性回归	$2.519\times10^{-4}$	12.24%
GM(1,1)	$5.440\times10^{-4}$	28.42%
Holt 平滑	$2.263\times10^{-4}$	10.77%
ARIMA	$2.951\times10^{-4}$	16.09%

图 2 以柱状展示四法 RMSE，图 7、图 8 进一步给出残差分布。综合看，Holt 指数平滑在留出检验中精度最高（MAPE 最低，约 10.8%），线性回归紧随其后（12.2%）；而 GM(1,1) 因指数结构在外推两步时放大误差，MAPE 高达 28.4%，表现最弱。这说明在"趋势明确但含噪声、样本少"的古塔变形场景下，稳健捕捉线性趋势的方法（Holt、线性）优于强假设的指数外推（GM），也优于依赖差分平稳性的 ARIMA。需要补充的是，留出验证采用"训练 8 期、预测后 2 期"的设定，比"全量训练预测未来 3 期"更严苛，也更贴近工程实际：决策正是要用已知历史预测未知近期，而非用全部已知回测。因此本文以留出 RMSE/MAPE 为方法选型的硬指标，样本内拟合优度仅作参考——后者易因过拟合而误导。

图 7 给出 t=8、t=9 两个留出点上各方法的绝对残差：Holt 平滑在两点均保持最小量级，GM(1,1) 在 t=9 处残差显著放大，印证了指数外推在多步预测下误差累积的机理。图 8 以预测值—实际值散点作对角线检验，点越贴近 45° 线说明预测越准，Holt 与线性回归的点最贴合，GM(1,1) 明显偏离，与 RMSE/MAPE 排序完全一致，构成对方法选型的第二重独立佐证。

八、未来三期预测对比

图 6 汇总四法对 t=10,11,12 的预测：线性回归、Holt 高度一致（约 0.00206→0.00228），GM(1,1) 略高（0.00216→0.00252），ARIMA 最低（约 0.00177）。以 ARIMA（保守低估）与 GM(1,1)（激进高估）作为预测区间的上下界参考，差距随期次拉大，反映外推越远不确定性越高。工程上取线性/Holt 作为点估计最稳妥。从区间宽度看，t=10 时四法最大最小差约 $4\times10^{-4}$ ，t=12 时扩大至约 $7.6\times10^{-4}$ ，说明外推越远、方法分歧越大，不确定性随之上升。这提示监测应"滚动预测、定期更新"：每获得新一期实测值就重新训练并外推，而非一次性外推多年，以免误差在远端累积放大。

九、方法选型讨论

当数据呈清晰单调趋势且样本稀少时，线性回归与 Holt 指数平滑是首选，二者稳健且可解释；GM(1,1) 虽能建模指数趋势，但本场景留出处误差最大，仅在确有指数加速且不做多步外推时谨慎使用。
Holt 指数平滑在保留趋势的同时自适应近期权重，与线性结果互证，是留出验证中的最优方法。
ARIMA 更适合平稳或周期序列；对强趋势短序列，差分会削弱趋势信号，需升阶（如 ARIMA(1,2,0)）或先去趋势再建模。
留出验证是方法选型的"试金石"：必须以对未知期的预测误差而非样本内拟合优度为准。
无论选何种方法，都应坚持"多模型互证"：用一种方法做点估计、其余方法做区间包络，并对预测结果标注"基于历史趋势外推、未考虑突发加固或灾害"的前提。预测本质是趋势延伸而非物理因果推断，这一点在给文保部门出具报告时必须明示。

十、结论

以倾斜量为对象的四种预测方法均给出"变形持续增大"的一致结论；留出验证表明 Holt 指数平滑精度最优（MAPE 10.8%）、线性回归稳健（12.2%），GM(1,1) 因指数外推过激进误差最大（28.4%）。点估计取线性与 Holt 一致值 $\hat{\tan\theta}_{10}\approx0.00206$ ，对应顶层水平偏移约 $0.20$ m（以 96 m 高计），已较 t=0（0.10 m）翻倍，警示未来需重点监测。范文三将把预测与四项指标的综合评价、灵敏度及蒙特卡洛区间结合，建立预警体系。综合而言，古塔倾斜预测不存在"唯一最优"的通用模型，只有"最适配数据形态"的模型：趋势平稳用线性/Holt，疑似指数加速可试 GM，平稳周期用 ARIMA。本文在给定合成数据下以 Holt 与线性回归为推荐使用，并以四法包络给出预测区间，形成可解释、可复核、可更新的预测流程。

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

import csv, math

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(p[0] 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 out

def gm11(x):
    n = len(x)
    x1 = [sum(x[:k + 1]) for k in range(n)]
    B = [[-0.5 * (x1[k] + x1[k - 1]), 1.0] for k in range(1, n)]
    Y = [x[k] for k in range(1, n)]
    m00 = sum(B[i][0] * B[i][0] for i in range(len(B)))
    m01 = sum(B[i][0] * B[i][1] for i in range(len(B)))
    m11 = sum(B[i][1] * B[i][1] for i in range(len(B)))
    y0 = sum(B[i][0] * Y[i] for i in range(len(B)))
    y1 = sum(B[i][1] * Y[i] for i in range(len(B)))
    det = m00 * m11 - m01 * m01
    a = (m11 * y0 - m01 * y1) / det
    u = (m00 * y1 - m01 * y0) / det
    x10 = x[0] - u / a
    def fc(h):
        return [x10 * math.exp(-a * k) + u / a - (x10 * math.exp(-a * (k - 1)) + u / a)
                for k in range(n, n + h)]
    return fc

def holt(x, alpha=0.3, beta=0.1):
    n = len(x); level = x[0]; trend = x[1] - x[0] if n > 1 else 0
    fit = [level]
    for i in range(1, n):
        last = level
        level = alpha * x[i] + (1 - alpha) * (level + trend)
        trend = beta * (level - last) + (1 - beta) * trend
        fit.append(level)
    def fc(h):
        return [level + trend * (k + 1) for k in range(h)]
    return fc

def arima1110(x):
    d = [x[i] - x[i - 1] for i in range(1, len(x))]
    num = sum(d[i] * d[i - 1] for i in range(1, len(d)))
    den = sum(d[i - 1] ** 2 for i in range(1, len(d)))
    phi = num / den if den else 0.0
    def fc(h):
        out = []; prev = x[-1]; dd = d[-1]
        for _ in range(h):
            dd = phi * dd; prev = prev + dd; out.append(prev)
        return out
    return fc

def rmse(a, b):
    return math.sqrt(sum((p - q) ** 2 for p, q in zip(a, b)) / len(a))
def mape(a, b):
    return sum(abs(p - q) / abs(q) for p, q in zip(a, b)) / len(a) * 100

data = load("cumcm2013c.csv")
C = centers(data)
tilt = tilt_series(C)
xs = list(range(len(tilt)))

a, b = linreg(xs, tilt)
fc_lin = [a * x + b for x in range(len(tilt), len(tilt) + 3)]
fc_g = gm11(tilt)(3)
fc_h = holt(tilt)(3)
fc_a = arima1110(tilt)(3)
print("线性回归 a,b :", round(a, 8), round(b, 8))
print("线性外推     :", [round(v, 6) for v in fc_lin])
print("GM(1,1)     :", [round(v, 6) for v in fc_g])
print("Holt        :", [round(v, 6) for v in fc_h])
print("ARIMA       :", [round(v, 6) for v in fc_a])

trx, tey = xs[:8], tilt[:8]
a2, b2 = linreg(trx, tey)
pv_lin = [a2 * x + b2 for x in [8, 9]]
pv_g = gm11(tey)(2)
pv_h = holt(tey)(2)
pv_a = arima1110(tey)(2)
actual = tilt[8:10]
for name, pv in [("线性回归", pv_lin), ("GM(1,1)", pv_g), ("Holt", pv_h), ("ARIMA", pv_a)]:
    print("%-8s RMSE=%.3e MAPE=%.3f%%" % (name, rmse(pv, actual), mape(pv, actual)))