结构健康监测仿真 - 主题037：结构健康监测中的不确定性量化技术

目录

1. 引言
2. 不确定性来源与分类
3. 概率不确定性量化方法
4. 非概率不确定性量化方法
5. 混合不确定性量化方法
6. 不确定性在SHM中的应用
7. Python仿真实践
8. 案例研究
9. 总结与展望

---

引言

1.1 不确定性量化的重要性

结构健康监测（Structural Health Monitoring, SHM）系统的核心任务是通过传感器数据评估结构的安全状态。然而，在实际工程应用中，从数据采集到损伤评估的每个环节都充满了各种不确定性。这些不确定性如果不加以量化和处理，可能导致错误的诊断决策，甚至引发严重的安全事故。

不确定性量化的重要性体现在：

• 风险评估：提供可靠的安全裕度评估
• 决策支持：为维修和维护决策提供置信度信息
• 模型验证：评估数值模型与真实结构的吻合程度
• 传感器优化：指导传感器布置和数据采集策略
• 标准制定：为结构安全评估标准提供科学依据

1.2 结构健康监测中的不确定性挑战

结构健康监测面临的不确定性挑战包括：

（1）认知不确定性（Epistemic Uncertainty）

• 模型简化导致的误差
• 参数识别的不完整性
• 边界条件的未知性
• 损伤机理的认知局限

（2）偶然不确定性（Aleatory Uncertainty）

• 材料属性的固有变异性
• 环境荷载的随机性
• 测量噪声的随机性
• 制造和施工偏差

（3）数据不确定性

• 传感器精度限制
• 采样频率不足
• 数据缺失和异常
• 传输和存储误差

（4）模型不确定性

• 物理模型简化
• 数值离散化误差
• 本构关系不准确
• 边界条件简化

1.3 不确定性量化的价值

通过系统的不确定性量化，可以实现：

• 提高诊断可靠性：给出损伤检测的置信度
• 优化资源配置：基于不确定性评估分配监测资源
• 延长结构寿命：在不确定性约束下优化维护策略
• 降低维护成本：避免过度维修和维修不足
• 提升安全水平：量化安全裕度，预防灾难性失效

---

不确定性来源与分类

2.1 不确定性的数学定义

随机变量：设$X$是定义在概率空间$(\Omega ,\mathcal{F},P)$上的函数，如果对任意实数$x$，事件 { ω : X ( ω ) ≤ x } ∈ F \{\omega: X(\omega) \leq x\} \in \mathcal{F} {ω:X(ω)≤x}∈F，则称$X$为随机变量。

不确定性量化：给定模型$Y=f(X)$，其中$X$为输入不确定性，量化输出$Y$的不确定性分布特征。

2.2 不确定性的分类体系

（1）按来源分类

类型	描述	示例
参数不确定性	模型参数的变异性	材料弹性模量、密度
模型不确定性	模型形式的误差	简化假设、本构关系
测量不确定性	数据采集的误差	传感器噪声、校准误差
边界不确定性	边界条件的未知性	支撑条件、连接刚度

（2）按性质分类

偶然不确定性（Aleatory）：

• 本质随机，不可消除
• 可用概率分布描述
• 例：风速、地震动、材料强度

认知不确定性（Epistemic）：

• 源于知识不足，可减少
• 可用区间、证据理论描述
• 例：模型误差、参数识别误差

（3）按传播方式分类

• 前向不确定性传播：输入→输出的不确定性传递
• 逆向不确定性量化：输出→输入的不确定性反演
• 模型选择不确定性：多个竞争模型的不确定性

2.3 不确定性的数学表征

（1）概率表征

概率密度函数（PDF）：$p_X(x)$

累积分布函数（CDF）：$F_X(x)=P(X\le x)={\int }_{-\infty }^x{p}_X(t)dt$

统计矩：

• 均值：${\mu }_X=E[X]={\int }_{-\infty }^{\infty }x{p}_X(x)dx$
• 方差：${\sigma }_X^2=E[(X-{\mu }_X{)}^2]$
• 偏度：${\gamma }_1=E[(X-{\mu }_X{)}^3]/{\sigma }_X^3$
• 峰度：${\gamma }_2=E[(X-{\mu }_X{)}^4]/{\sigma }_X^4-3$

（2）区间表征

区间变量：$X^I=[\underline{x},\overline{x}]$

区间运算：

• 加法：$[a,b]+[c,d]=[a+c,b+d]$
• 乘法：$[a,b]\times [c,d]=[\min (ac,ad,bc,bd),\max (ac,ad,bc,bd)]$

（3）模糊数表征

模糊集合：$\tilde{A}=\{(x,{\mu }_A(x))|x\in X\}$

隶属函数：${\mu }_A:X\to [0,1]$

三角模糊数：$\tilde{A}=(a,b,c)$，其中$b$为核心，$a,c$为边界

2.4 不确定性传播理论

（1）泰勒展开法

一阶近似：$Y\approx f({\mu }_X)+{\sum }_{i=1}^n\frac{\partial f}{\partial X_i}{|}_{{\mu }_X}(X_i-{\mu }_{X_i})$

均值：${\mu }_Y\approx f({\mu }_X)$

方差：${\sigma }_Y^2\approx {\sum }_{i=1}^n{\sum }_{j=1}^n\frac{\partial f}{\partial X_i}\frac{\partial f}{\partial X_j}Cov(X_i,X_j)$

（2）蒙特卡洛模拟

基本步骤：

1. 从输入分布$p_X(x)$采样$N$个样本 { x ( i ) } i = 1 N \{x^{(i)}\}_{i=1}^{N} {x(i)}i=1N​
2. 计算输出样本 { y ( i ) = f ( x ( i ) ) } i = 1 N \{y^{(i)} = f(x^{(i)})\}_{i=1}^{N} {y(i)=f(x(i))}i=1N​
3. 统计分析输出样本

收敛性：$\text{Var}({\hat{\mu }}_Y)={\sigma }_Y^2/N$，误差$O(1/\sqrt{N})$

（3）多项式混沌展开

将随机输出展开为随机变量的多项式：

$$Y(\xi )=\sum _{\alpha \in \mathcal{J}}{y}_{\alpha }{\Psi }_{\alpha }(\xi )$$

其中${\Psi }_{\alpha }$为多维正交多项式，$\xi$为标准随机变量。

---

概率不确定性量化方法

3.1 蒙特卡洛方法

基本原理：

利用大数定律，通过大量随机采样估计统计量：

$${\hat{\mu }}_Y=\frac{1}{N}\sum _{i=1}^N{Y}^{(i)}$$

方差缩减技术：

（1）重要性采样

改变采样分布，增加重要区域的采样密度：

$$E[g(X)]=\int g(x)\frac{p(x)}{q(x)}q(x)dx$$

（2）分层采样

将输入空间划分为若干层，每层独立采样：

$${\hat{\mu }}_Y=\sum _{j=1}^L{w}_j{\hat{\mu }}_j$$

（3）拉丁超立方采样（LHS）

确保每个维度都被均匀采样：

1. 将每个维度划分为$N$个等概率区间
2. 在每个维度上随机选择区间
3. 组合形成样本点

优点：

• 收敛速度与维度无关
• 实现简单，适用性广
• 可并行计算

缺点：

• 收敛速度慢（$O(1/\sqrt{N})$）
• 高维问题计算成本高
• 小概率事件估计困难

3.2 多项式混沌展开（PCE）

基本理论：

将随机输出展开为随机输入的正交多项式：

$$Y(\xi )=\sum _{\alpha \in \mathcal{A}}{c}_{\alpha }{\Psi }_{\alpha }(\xi )$$

其中${\Psi }_{\alpha }$为关于概率测度的正交多项式。

常用正交多项式：

分布类型	多项式类型	权重函数
高斯分布	Hermite	$e^{-x^2/2}$
均匀分布	Legendre	1
Gamma分布	Laguerre	$x^{\alpha }{e}^{-x}$
Beta分布	Jacobi	$(1-x{)}^{\alpha }(1+x{)}^{\beta }$

系数计算方法：

（1）投影法

$$c_{\alpha }=\frac{E[Y{\Psi }_{\alpha }]}{E[{\Psi }_{\alpha }^2]}=\frac{1}{{\gamma }_{\alpha }}\int Y(\xi ){\Psi }_{\alpha }(\xi )p(\xi )d\xi$$

使用高斯积分或蒙特卡洛积分计算。

（2）回归法

通过最小二乘拟合：

min ⁡ { c α } ∑ i = 1 N ( Y ( i ) − ∑ α c α Ψ α ( ξ ( i ) ) ) 2 \min_{\{c_\alpha\}} \sum_{i=1}^{N} \left(Y^{(i)} - \sum_{\alpha} c_\alpha \Psi_\alpha(\xi^{(i)})\right)^2 {cα​}min​i=1∑N​(Y(i)−α∑​cα​Ψα​(ξ(i)))2

统计量计算：

均值：${\mu }_Y=c_0$

方差：${\sigma }_Y^2={\sum }_{\alpha \ne 0}{c}_{\alpha }^2{\gamma }_{\alpha }$

优点：

• 收敛速度快（指数收敛）
• 计算效率高（一次构建，多次使用）
• 提供敏感性分析

缺点：

• 高维问题存在"维度灾难"
• 非光滑问题收敛慢
• 需要知道输入分布

3.3 随机配点法

基本思想：

在精心选择的配点上求解确定性问题，然后插值得到统计量。

配点选择：

（1）张量积配点

一维配点的笛卡尔积： { ξ i } × { ξ j } \{\xi_i\} \times \{\xi_j\} {ξi​}×{ξj​}

（2）稀疏网格配点

Smolyak算法减少配点数量：

$$A(q,n)=\sum _{q-n+1\le |\mathbf{i}|\le q}(-1{)}^{q-|\mathbf{i}|}(\genfrac{}{}{0ex}{}{n-1}{q-|\mathbf{i}|})(U^{i_1}\otimes \cdots \otimes U^{i_n})$$

插值方法：

拉格朗日插值：$Y(\xi )\approx {\sum }_{i=1}^{N}Y({\xi }^{(i)})L_i(\xi )$

优点：

• 计算效率高于蒙特卡洛
• 适合光滑问题
• 可直接使用确定性求解器

缺点：

• 高维问题配点数量指数增长
• 非光滑问题精度低
• 需要自适应策略

3.4 随机有限元方法

摄动法：

对随机输入进行泰勒展开：

$$K(\theta )=K_0+\sum _i\frac{\partial K}{\partial {\theta }_i}\Delta {\theta }_i+\frac{1}{2}\sum _{i,j}\frac{{\partial }^{2}K}{\partial {\theta }_i\partial {\theta }_j}\Delta {\theta }_i\Delta {\theta }_j+\cdots$$

Neumann展开：

$$(K_0+\Delta K)u=f$$

$$u=(I+K_0^{-1}\Delta K{)}^{-1}{K}_0^{-1}f\approx \sum _{i=0}^m(-K_0^{-1}\Delta K{)}^i{u}_0$$

谱随机有限元：

将随机场离散为随机变量：

$$E(x,\theta )=\bar{E}(x)+\sum _{i=1}^M\sqrt{{\lambda }_i}{\varphi }_i(x){\xi }_i(\theta )$$

其中$({\lambda }_i,{\varphi }_i)$为协方差核的特征对。

---

非概率不确定性量化方法

4.1 区间分析方法

基本定义：

区间数： [ a ] = [ a ‾ , a ‾ ] = { x ∈ R ∣ a ‾ ≤ x ≤ a ‾ } [a] = [\underline{a}, \overline{a}] = \{x \in \mathbb{R} | \underline{a} \leq x \leq \overline{a}\} [a]=[a​,a]={x∈R∣a​≤x≤a}

区间宽度：$w([a])=\overline{a}-\underline{a}$

区间中点：$m([a])=(\underline{a}+\overline{a})/2$

区间运算：

加法：$[a]+[b]=[\underline{a}+\underline{b},\overline{a}+\overline{b}]$

减法：$[a]-[b]=[\underline{a}-\overline{b},\overline{a}-\underline{b}]$

乘法：$[a]\times [b]=[\min (S),\max (S)]$，其中 S = { a ‾ b ‾ , a ‾ b ‾ , a ‾ b ‾ , a ‾ b ‾ } S = \{\underline{a}\underline{b}, \underline{a}\overline{b}, \overline{a}\underline{b}, \overline{a}\overline{b}\} S={a​b​,a​b,ab​,ab}

除法：$[a]/[b]=[a]\times [1/\overline{b},1/\underline{b}]$，要求$0\notin [b]$

区间有限元：

$$[K][u]=[f]$$

求解区间线性方程组，得到响应区间。

仿射算术：

减少区间扩张：$\hat{x}=x_0+{\sum }_{i=1}^n{x}_i{\epsilon }_i$，其中${\epsilon }_i\in [-1,1]$

4.2 证据理论（Dempster-Shafer理论）

基本框架：

识别框架： Θ = { θ 1 , θ 2 , . . . , θ n } \Theta = \{\theta_1, \theta_2, ..., \theta_n\} Θ={θ1​,θ2​,...,θn​}

基本概率分配（BPA）：$m:2^{\Theta }\to [0,1]$，满足$m(\varnothing )=0$，${\sum }_{A\subseteq \Theta }m(A)=1$

信任函数与似然函数：

信任函数：$Bel(A)={\sum }_{B\subseteq A}m(B)$

似然函数：$Pl(A)={\sum }_{B\cap A\ne \varnothing }m(B)$

Dempster组合规则：

$$m_{1,2}(A)=\frac{1}{1-K}\sum _{B\cap C=A}{m}_1(B)m_2(C)$$

其中$K={\sum }_{B\cap C=\varnothing }{m}_1(B)m_2(C)$

在SHM中的应用：

• 多源证据融合
• 专家知识整合
• 不确定性传播

4.3 可能性理论

模糊集合：

隶属函数：${\mu }_A:X\to [0,1]$

可能性分布：

$\pi (x)={\mu }_A(x)$，表示$x$是真实值的可能性

可能性测度：

$\Pi (A)={\sup }_{x\in A}\pi (x)$

必要性测度：

$N(A)=1-\Pi (A^c)={\inf }_{x\notin A}(1-\pi (x))$

模糊运算：

扩展原理：${\mu }_{f(A)}(y)={\sup }_{x:f(x)=y}{\mu }_A(x)$

在SHM中的应用：

• 语言变量建模
• 专家经验量化
• 模糊损伤识别

4.4 概率边界分析（P-Box）

定义：

概率边界：$[\underline{F}(x),\overline{F}(x)]$，其中$\underline{F}(x)\le F(x)\le \overline{F}(x)$

构建方法：

（1）置信区间法

基于Kolmogorov-Smirnov统计量：

$${\underline{F}}_n(x)=\max \left(0,F_n(x)-D_{n,\alpha }\right)$$

$${\overline{F}}_n(x)=\min \left(1,F_n(x)+D_{n,\alpha }\right)$$

（2）分布参数不确定

当分布类型已知但参数不确定时：

$$\underline{F}(x)=\underset{\theta \in [\underline{\theta },\overline{\theta }]}{\inf }F(x;\theta )$$

$$\overline{F}(x)=\underset{\theta \in [\underline{\theta },\overline{\theta }]}{\sup }F(x;\theta )$$

传播方法：

区间蒙特卡洛：在每个样本点计算输出边界

---

混合不确定性量化方法

5.1 概率-区间混合方法

问题描述：

输入$X=(X_p,X_i)$，其中$X_p$为概率变量，$X_i$为区间变量。

双层蒙特卡洛：

外层：对$X_p$采样
内层：对$X_i$进行区间分析

$$\underline{Y}(X_p)=\underset{X_i}{\min }f(X_p,X_i)$$

$$\overline{Y}(X_p)=\underset{X_i}{\max }f(X_p,X_i)$$

效率优化：

使用代理模型替代内层优化：

$$\tilde{Y}=\text{Surrogate}(X_p,X_i)$$

5.2 概率-证据混合方法

框架：

概率变量描述偶然不确定性
证据变量描述认知不确定性

传播策略：

1. 对概率变量进行蒙特卡洛采样
2. 对每个样本进行证据推理
3. 聚合所有样本的结果

应用：

• 模型参数（概率）+ 模型选择（证据）
• 材料属性（概率）+ 边界条件（证据）

5.3 随机区间有限元

基本方程：

$$K(\theta ,[\eta ])u=f(\theta )$$

其中$\theta$为随机变量，$[\eta ]$为区间变量。

求解策略：

（1）顺序求解

先进行随机分析，再进行区间分析

（2）联合展开

将随机和区间变量统一展开：

$$u(\theta ,\eta )=\sum _{i,j}{c}_{ij}{\Psi }_i(\theta ){\varphi }_j(\eta )$$

5.4 贝叶斯模型选择

模型选择不确定性：

当存在多个竞争模型$M_1,M_2,...,M_k$时，使用贝叶斯方法量化模型概率。

贝叶斯模型平均：

$$p(Y|D)=\sum _{i=1}^{k}p(Y|M_i,D)p(M_i|D)$$

模型证据：

$$p(D|M_i)=\int p(D|{\theta }_i,M_i)p({\theta }_i|M_i)d{\theta }_i$$

模型后验概率：

$$p(M_i|D)=\frac{p(D|M_i)p(M_i)}{{\sum }_{j}p(D|M_j)p(M_j)}$$

---

不确定性在SHM中的应用

6.1 损伤检测的不确定性量化

问题描述：

给定特征$D$，判断结构状态 H ∈ { H 0 , H 1 } H \in \{H_0, H_1\} H∈{H0​,H1​}（无损伤/有损伤）

贝叶斯方法：

$$P(H_1|D)=\frac{p(D|H_1)P(H_1)}{p(D|H_0)P(H_0)+p(D|H_1)P(H_1)}$$

决策规则：

若$P(H_1|D)>\eta$，则判定为损伤，其中$\eta$为决策阈值。

错误率估计：

虚警率：$P_{FA}=P(\text{判定损伤}|H_0)$

漏检率：$P_{MD}=P(\text{判定无损伤}|H_1)$

6.2 模型修正的不确定性

问题描述：

修正模型参数$\theta$使得模型响应$Y_m$与实测响应$Y_e$一致。

贝叶斯模型修正：

$$p(\theta |Y_e)\propto p(Y_e|\theta )p(\theta )$$

不确定性来源：

• 测量噪声
• 模型误差
• 参数相关性
• 可识别性问题

后验分析：

• 参数边缘分布
• 置信区间
• 相关系数矩阵

6.3 剩余寿命预测的不确定性

退化模型：

$$D(t)=D_0+{\int }_0^{t}r(s)ds+\epsilon (t)$$

其中$r(t)$为退化速率，$\epsilon (t)$为随机噪声。

不确定性来源：

• 初始状态不确定性
• 退化速率变异性
• 环境荷载随机性
• 模型误差

预测区间：

使用蒙特卡洛或多项式混沌预测失效时间分布：

$$P(T_f\le t)=P(D(t)\ge D_{threshold})$$

6.4 传感器优化布置

信息熵准则：

最大化传感器布置的信息增益：

$$\underset{S}{\max }H(Y)-H(Y|S)$$

其中$S$为传感器布置方案。

不确定性考虑：

• 测量噪声水平
• 传感器失效概率
• 环境干扰
• 成本约束

---

Python仿真实践

7.1 仿真目标

本仿真实现一个综合的不确定性量化系统，包括：

• 蒙特卡洛模拟
• 多项式混沌展开
• 区间分析
• 贝叶斯推断
• 不确定性传播可视化

7.2 完整Python代码

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
import imageio
import os
from scipy import stats
from scipy.integrate import quad
from scipy.interpolate import interp1d
from scipy.optimize import minimize

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

def generate_structural_response(E_mean=200e9, E_std=10e9, rho_mean=7850, rho_std=200, 
                                   L=10.0, A=0.01, F=1000, n_samples=1000):
    """生成结构响应（悬臂梁端部位移）
    
    位移公式：δ = FL³/(3EI)，其中I = Ah²/12（矩形截面）
    假设h = sqrt(A)，则 I = A²/12
    """
    np.random.seed(42)
    
    # 随机采样
    E_samples = np.random.normal(E_mean, E_std, n_samples)
    rho_samples = np.random.normal(rho_mean, rho_std, n_samples)
    
    # 计算截面惯性矩（假设正方形截面）
    h = np.sqrt(A)  # 截面高度
    I = A * h**2 / 12  # 惯性矩
    
    # 计算位移（考虑材料不确定性）
    displacement = F * L**3 / (3 * E_samples * I)
    
    # 添加模型误差（认知不确定性）
    model_error = np.random.uniform(0.95, 1.05, n_samples)
    displacement_with_error = displacement * model_error
    
    return displacement_with_error, E_samples, rho_samples

def monte_carlo_simulation(n_samples=10000):
    """蒙特卡洛模拟"""
    print('\n执行蒙特卡洛模拟...')
    
    displacement, E_samples, rho_samples = generate_structural_response(n_samples=n_samples)
    
    # 统计量
    mean_disp = np.mean(displacement)
    std_disp = np.std(displacement)
    ci_lower = np.percentile(displacement, 2.5)
    ci_upper = np.percentile(displacement, 97.5)
    
    print(f'   样本数: {n_samples}')
    print(f'   位移均值: {mean_disp:.6f} m')
    print(f'   位移标准差: {std_disp:.6f} m')
    print(f'   95%置信区间: [{ci_lower:.6f}, {ci_upper:.6f}] m')
    
    return displacement, mean_disp, std_disp, (ci_lower, ci_upper)

def polynomial_chaos_expansion(order=3, n_samples=1000):
    """多项式混沌展开（简化版，使用Legendre多项式）"""
    print('\n执行多项式混沌展开...')
    
    # 定义Legendre多项式（简化，使用numpy的legendre）
    def legendre_poly(x, n):
        """计算n阶Legendre多项式"""
        if n == 0:
            return np.ones_like(x)
        elif n == 1:
            return x
        elif n == 2:
            return 0.5 * (3*x**2 - 1)
        elif n == 3:
            return 0.5 * (5*x**3 - 3*x)
        else:
            return np.ones_like(x)
    
    # 将弹性模量标准化到[-1, 1]
    E_mean, E_std = 200e9, 10e9
    E_normalized = lambda xi: E_mean + E_std * xi * np.sqrt(3)
    
    # 生成训练样本
    np.random.seed(42)
    xi_train = np.random.uniform(-1, 1, n_samples)
    E_train = E_normalized(xi_train)
    
    # 计算响应
    L, A, F = 10.0, 0.01, 1000
    h = np.sqrt(A)
    I = A * h**2 / 12
    displacement_train = F * L**3 / (3 * E_train * I)
    
    # 构建设计矩阵
    Phi = np.zeros((n_samples, order+1))
    for i in range(order+1):
        Phi[:, i] = legendre_poly(xi_train, i)
    
    # 最小二乘求解系数
    coeffs, residuals, rank, s = np.linalg.lstsq(Phi, displacement_train, rcond=None)
    
    print(f'   展开阶数: {order}')
    print(f'   PCE系数: {coeffs}')
    print(f'   均值估计: {coeffs[0]:.6f} m')
    print(f'   方差估计: {np.sum(coeffs[1:]**2 / (2*np.arange(1,order+1)+1)):.6e} m²')
    
    return coeffs, xi_train, displacement_train

def interval_analysis():
    """区间分析"""
    print('\n执行区间分析...')
    
    # 定义区间变量
    E_interval = [190e9, 210e9]  # 弹性模量区间
    F_interval = [900, 1100]      # 荷载区间
    
    L, A = 10.0, 0.01
    h = np.sqrt(A)
    I = A * h**2 / 12
    
    # 计算位移区间（使用区间运算）
    # δ = FL³/(3EI)
    # 对于单调函数，极值出现在边界点
    
    combinations = [
        (E_interval[0], F_interval[0]),
        (E_interval[0], F_interval[1]),
        (E_interval[1], F_interval[0]),
        (E_interval[1], F_interval[1])
    ]
    
    displacements = []
    for E, F in combinations:
        disp = F * L**3 / (3 * E * I)
        displacements.append(disp)
    
    disp_lower = min(displacements)
    disp_upper = max(displacements)
    
    print(f'   弹性模量区间: [{E_interval[0]/1e9:.1f}, {E_interval[1]/1e9:.1f}] GPa')
    print(f'   荷载区间: [{F_interval[0]}, {F_interval[1]}] N')
    print(f'   位移区间: [{disp_lower:.6f}, {disp_upper:.6f}] m')
    print(f'   区间宽度: {disp_upper - disp_lower:.6f} m')
    
    return (disp_lower, disp_upper), E_interval, F_interval

def bayesian_inference(n_observations=50):
    """贝叶斯推断（估计弹性模量）"""
    print('\n执行贝叶斯推断...')
    
    # 真实参数
    E_true = 200e9
    L, A, F = 10.0, 0.01, 1000
    h = np.sqrt(A)
    I = A * h**2 / 12
    
    # 生成观测数据（带噪声）
    np.random.seed(42)
    disp_true = F * L**3 / (3 * E_true * I)
    noise_std = 0.0001  # 测量噪声标准差
    observations = disp_true + np.random.normal(0, noise_std, n_observations)
    
    # 定义先验（均匀分布）
    E_prior_range = [180e9, 220e9]
    
    # 定义似然函数
    def likelihood(E, obs):
        disp_model = F * L**3 / (3 * E * I)
        return np.prod(stats.norm.pdf(obs, disp_model, noise_std))
    
    # 网格搜索计算后验
    E_grid = np.linspace(E_prior_range[0], E_prior_range[1], 1000)
    posteriors = []
    
    for E in E_grid:
        post = likelihood(E, observations)
        posteriors.append(post)
    
    posteriors = np.array(posteriors)
    posteriors = posteriors / np.sum(posteriors)  # 归一化
    
    # 计算后验统计量
    E_map = E_grid[np.argmax(posteriors)]
    E_mean_post = np.sum(E_grid * posteriors)
    E_std_post = np.sqrt(np.sum((E_grid - E_mean_post)**2 * posteriors))
    
    # 95%可信区间
    cdf = np.cumsum(posteriors)
    ci_lower = E_grid[np.searchsorted(cdf, 0.025)]
    ci_upper = E_grid[np.searchsorted(cdf, 0.975)]
    
    print(f'   观测次数: {n_observations}')
    print(f'   MAP估计: {E_map/1e9:.2f} GPa')
    print(f'   后验均值: {E_mean_post/1e9:.2f} GPa')
    print(f'   后验标准差: {E_std_post/1e9:.2f} GPa')
    print(f'   95%可信区间: [{ci_lower/1e9:.2f}, {ci_upper/1e9:.2f}] GPa')
    
    return E_grid, posteriors, observations, (E_map, E_mean_post, E_std_post)

def sensitivity_analysis():
    """敏感性分析（Sobol指数）"""
    print('\n执行敏感性分析...')
    
    # 定义参数范围
    E_mean, E_std = 200e9, 10e9
    rho_mean, rho_std = 7850, 200
    F_mean, F_std = 1000, 50
    
    L, A = 10.0, 0.01
    h = np.sqrt(A)
    I = A * h**2 / 12
    
    # 蒙特卡洛估计Sobol指数（简化版）
    n_samples = 5000
    np.random.seed(42)
    
    # 采样
    E_samples = np.random.normal(E_mean, E_std, n_samples)
    rho_samples = np.random.normal(rho_mean, rho_std, n_samples)
    F_samples = np.random.normal(F_mean, F_std, n_samples)
    
    # 计算响应
    displacement = F_samples * L**3 / (3 * E_samples * I)
    
    # 计算总方差
    total_var = np.var(displacement)
    
    # 一阶效应（主效应）
    # 固定E，变化其他参数
    E_fixed = E_mean
    disp_E_fixed = F_samples * L**3 / (3 * E_fixed * I)
    var_without_E = np.var(disp_E_fixed)
    S_E = (total_var - var_without_E) / total_var
    
    # 固定F，变化其他参数
    F_fixed = F_mean
    disp_F_fixed = F_fixed * L**3 / (3 * E_samples * I)
    var_without_F = np.var(disp_F_fixed)
    S_F = (total_var - var_without_F) / total_var
    
    print(f'   总方差: {total_var:.6e}')
    print(f'   弹性模量Sobol指数: {S_E:.4f}')
    print(f'   荷载Sobol指数: {S_F:.4f}')
    print(f'   密度Sobol指数: {0:.4f} (对位移无直接影响)')
    
    return {'E': S_E, 'F': S_F, 'rho': 0.0}, total_var

def reliability_analysis():
    """可靠性分析"""
    print('\n执行可靠性分析...')
    
    # 定义极限状态函数
    # g(X) = δ_allowable - δ_actual
    # 失效当g(X) < 0
    
    delta_allowable = 0.0009  # 允许位移
    
    E_mean, E_std = 200e9, 10e9
    F_mean, F_std = 1000, 50
    L, A = 10.0, 0.01
    h = np.sqrt(A)
    I = A * h**2 / 12
    
    # 蒙特卡洛可靠性分析
    n_samples = 100000
    np.random.seed(42)
    
    E_samples = np.random.normal(E_mean, E_std, n_samples)
    F_samples = np.random.normal(F_mean, F_std, n_samples)
    
    displacement = F_samples * L**3 / (3 * E_samples * I)
    g_values = delta_allowable - displacement
    
    # 失效概率
    failure_count = np.sum(g_values < 0)
    Pf = failure_count / n_samples
    
    # 可靠度指标
    beta = -stats.norm.ppf(Pf)
    
    print(f'   允许位移: {delta_allowable:.6f} m')
    print(f'   失效概率 Pf: {Pf:.6f}')
    print(f'   可靠度指标 β: {beta:.4f}')
    print(f'   失效次数: {failure_count}/{n_samples}')
    
    return Pf, beta, g_values

def visualize_uncertainty_analysis(save_path='uncertainty_analysis.png'):
    """可视化不确定性分析结果"""
    print('\n生成不确定性分析可视化...')
    
    fig = plt.figure(figsize=(20, 16))
    
    # 1. 蒙特卡洛结果直方图
    ax1 = plt.subplot(4, 3, 1)
    displacement, mean_disp, std_disp, ci = monte_carlo_simulation(n_samples=10000)
    ax1.hist(displacement, bins=50, density=True, alpha=0.7, color='skyblue', edgecolor='black')
    ax1.axvline(mean_disp, color='red', linestyle='--', linewidth=2, label=f'均值={mean_disp:.6f}')
    ax1.axvline(ci[0], color='orange', linestyle=':', linewidth=2, label=f'95% CI')
    ax1.axvline(ci[1], color='orange', linestyle=':', linewidth=2)
    ax1.set_xlabel('位移 (m)')
    ax1.set_ylabel('概率密度')
    ax1.set_title('蒙特卡洛模拟结果')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 位移累积分布函数
    ax2 = plt.subplot(4, 3, 2)
    sorted_disp = np.sort(displacement)
    cdf = np.arange(1, len(sorted_disp)+1) / len(sorted_disp)
    ax2.plot(sorted_disp, cdf, linewidth=2, color='blue')
    ax2.axhline(0.025, color='red', linestyle='--', alpha=0.5)
    ax2.axhline(0.975, color='red', linestyle='--', alpha=0.5)
    ax2.set_xlabel('位移 (m)')
    ax2.set_ylabel('累积概率')
    ax2.set_title('位移累积分布函数')
    ax2.grid(True, alpha=0.3)
    
    # 3. 多项式混沌展开
    ax3 = plt.subplot(4, 3, 3)
    coeffs, xi_train, disp_train = polynomial_chaos_expansion(order=3, n_samples=1000)
    xi_test = np.linspace(-1, 1, 100)
    
    # 计算PCE预测
    def legendre_poly(x, n):
        if n == 0:
            return np.ones_like(x)
        elif n == 1:
            return x
        elif n == 2:
            return 0.5 * (3*x**2 - 1)
        elif n == 3:
            return 0.5 * (5*x**3 - 3*x)
        else:
            return np.ones_like(x)
    
    disp_pce = np.zeros_like(xi_test)
    for i, c in enumerate(coeffs):
        disp_pce += c * legendre_poly(xi_test, i)
    
    ax3.scatter(xi_train[:100], disp_train[:100], alpha=0.5, s=20, label='样本')
    ax3.plot(xi_test, disp_pce, 'r-', linewidth=2, label='PCE近似')
    ax3.set_xlabel('标准化输入')
    ax3.set_ylabel('位移 (m)')
    ax3.set_title('多项式混沌展开')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 区间分析
    ax4 = plt.subplot(4, 3, 4)
    disp_interval, E_interval, F_interval = interval_analysis()
    
    # 绘制区间
    ax4.barh(['位移'], [disp_interval[1] - disp_interval[0]], 
             left=[disp_interval[0]], height=0.5, color='lightcoral', alpha=0.7, edgecolor='black')
    ax4.axvline(disp_interval[0], color='red', linestyle='--', linewidth=2)
    ax4.axvline(disp_interval[1], color='red', linestyle='--', linewidth=2)
    ax4.set_xlabel('位移 (m)')
    ax4.set_title('区间分析结果')
    ax4.grid(True, alpha=0.3, axis='x')
    
    # 5. 贝叶斯后验分布
    ax5 = plt.subplot(4, 3, 5)
    E_grid, posteriors, observations, stats = bayesian_inference(n_observations=50)
    ax5.plot(E_grid/1e9, posteriors, linewidth=2, color='blue')
    ax5.axvline(stats[0]/1e9, color='red', linestyle='--', linewidth=2, label=f'MAP={stats[0]/1e9:.1f} GPa')
    ax5.axvline(200, color='green', linestyle=':', linewidth=2, label='真实值=200 GPa')
    ax5.set_xlabel('弹性模量 (GPa)')
    ax5.set_ylabel('后验概率密度')
    ax5.set_title('贝叶斯后验分布')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 6. 敏感性分析
    ax6 = plt.subplot(4, 3, 6)
    sobol_indices, total_var = sensitivity_analysis()
    params = list(sobol_indices.keys())
    indices = list(sobol_indices.values())
    
    colors = ['skyblue', 'lightgreen', 'lightcoral']
    bars = ax6.bar(params, indices, color=colors, edgecolor='black')
    ax6.set_ylabel('Sobol指数')
    ax6.set_title('全局敏感性分析')
    ax6.set_ylim(0, 1)
    ax6.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar, idx in zip(bars, indices):
        height = bar.get_height()
        ax6.text(bar.get_x() + bar.get_width()/2., height,
                f'{idx:.3f}', ha='center', va='bottom')
    
    # 7. 可靠性分析
    ax7 = plt.subplot(4, 3, 7)
    Pf, beta, g_values = reliability_analysis()
    
    ax7.hist(g_values, bins=50, density=True, alpha=0.7, color='lightblue', edgecolor='black')
    ax7.axvline(0, color='red', linestyle='--', linewidth=2, label='极限状态')
    ax7.axvline(np.mean(g_values), color='green', linestyle='-', linewidth=2, 
                label=f'均值={np.mean(g_values):.6f}')
    ax7.set_xlabel('极限状态函数 g(X)')
    ax7.set_ylabel('概率密度')
    ax7.set_title(f'可靠性分析 (Pf={Pf:.6f}, β={beta:.2f})')
    ax7.legend()
    ax7.grid(True, alpha=0.3)
    
    # 8. 不确定性传播对比
    ax8 = plt.subplot(4, 3, 8)
    methods = ['蒙特卡洛', 'PCE', '区间分析']
    means = [mean_disp, coeffs[0], (disp_interval[0]+disp_interval[1])/2]
    stds = [std_disp, np.sqrt(np.sum(coeffs[1:]**2 / (2*np.arange(1,len(coeffs))+1))), 
            (disp_interval[1]-disp_interval[0])/4]
    
    x_pos = np.arange(len(methods))
    ax8.bar(x_pos, means, yerr=stds, capsize=5, color='skyblue', edgecolor='black', alpha=0.7)
    ax8.set_xticks(x_pos)
    ax8.set_xticklabels(methods)
    ax8.set_ylabel('位移 (m)')
    ax8.set_title('不同方法的不确定性量化对比')
    ax8.grid(True, alpha=0.3, axis='y')
    
    # 9. 参数相关性分析
    ax9 = plt.subplot(4, 3, 9)
    displacement, E_samples, rho_samples = generate_structural_response(n_samples=1000)
    
    # 计算相关系数
    corr_matrix = np.corrcoef([E_samples, rho_samples, displacement])
    im = ax9.imshow(corr_matrix, cmap='coolwarm', vmin=-1, vmax=1)
    ax9.set_xticks([0, 1, 2])
    ax9.set_yticks([0, 1, 2])
    ax9.set_xticklabels(['E', 'ρ', 'δ'])
    ax9.set_yticklabels(['E', 'ρ', 'δ'])
    ax9.set_title('参数相关性矩阵')
    plt.colorbar(im, ax=ax9)
    
    # 添加数值
    for i in range(3):
        for j in range(3):
            ax9.text(j, i, f'{corr_matrix[i,j]:.2f}', ha='center', va='center', 
                    color='white' if abs(corr_matrix[i,j]) > 0.5 else 'black')
    
    # 10. 收敛性分析
    ax10 = plt.subplot(4, 3, 10)
    sample_sizes = [100, 500, 1000, 5000, 10000]
    means_conv = []
    stds_conv = []
    
    for n in sample_sizes:
        disp, _, _ = generate_structural_response(n_samples=n)
        means_conv.append(np.mean(disp))
        stds_conv.append(np.std(disp))
    
    ax10.plot(sample_sizes, means_conv, 'o-', linewidth=2, markersize=8, label='均值')
    ax10.fill_between(sample_sizes, 
                       np.array(means_conv) - np.array(stds_conv),
                       np.array(means_conv) + np.array(stds_conv),
                       alpha=0.3, label='±1标准差')
    ax10.set_xlabel('样本数')
    ax10.set_ylabel('位移 (m)')
    ax10.set_title('蒙特卡洛收敛性')
    ax10.set_xscale('log')
    ax10.legend()
    ax10.grid(True, alpha=0.3)
    
    # 11. 认知vs偶然不确定性
    ax11 = plt.subplot(4, 3, 11)
    
    # 偶然不确定性（概率）
    aleatory_samples = np.random.normal(mean_disp, std_disp, 1000)
    
    # 认知不确定性（区间）
    epistemic_lower = mean_disp * 0.95
    epistemic_upper = mean_disp * 1.05
    
    ax11.hist(aleatory_samples, bins=30, density=True, alpha=0.5, 
              color='blue', label='偶然不确定性', edgecolor='black')
    ax11.axvspan(epistemic_lower, epistemic_upper, alpha=0.3, color='red', 
                 label='认知不确定性')
    ax11.axvline(mean_disp, color='black', linestyle='--', linewidth=2)
    ax11.set_xlabel('位移 (m)')
    ax11.set_ylabel('概率密度')
    ax11.set_title('认知vs偶然不确定性')
    ax11.legend()
    ax11.grid(True, alpha=0.3)
    
    # 12. 不确定性分解
    ax12 = plt.subplot(4, 3, 12)
    
    # 总不确定性分解
    categories = ['材料属性', '几何参数', '荷载', '模型误差', '测量误差']
    contributions = [0.4, 0.1, 0.25, 0.15, 0.1]  # 假设的贡献比例
    colors_pie = ['skyblue', 'lightgreen', 'lightcoral', 'gold', 'plum']
    
    wedges, texts, autotexts = ax12.pie(contributions, labels=categories, autopct='%1.1f%%',
                                         colors=colors_pie, startangle=90)
    ax12.set_title('不确定性来源分解')
    
    plt.tight_layout()
    plt.savefig(save_path, dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f'可视化已保存: {save_path}')

def create_uncertainty_animation(save_path='uncertainty_animation.gif'):
    """创建不确定性传播动画"""
    print('\n生成不确定性传播动画...')
    
    temp_files = []
    n_frames = 20
    
    # 参数范围
    E_mean, E_std = 200e9, 10e9
    L, A, F = 10.0, 0.01, 1000
    h = np.sqrt(A)
    I = A * h**2 / 12
    
    for frame_idx in range(n_frames):
        fig = plt.figure(figsize=(16, 10))
        
        # 逐渐增加样本数
        n_samples = (frame_idx + 1) * 100
        np.random.seed(42)
        
        E_samples = np.random.normal(E_mean, E_std, n_samples)
        displacement = F * L**3 / (3 * E_samples * I)
        
        # 1. 输入分布
        ax1 = fig.add_subplot(2, 3, 1)
        ax1.hist(E_samples/1e9, bins=20, density=True, alpha=0.7, color='skyblue', edgecolor='black')
        ax1.set_xlabel('弹性模量 (GPa)')
        ax1.set_ylabel('概率密度')
        ax1.set_title(f'输入分布 (n={n_samples})')
        ax1.grid(True, alpha=0.3)
        
        # 2. 输出分布
        ax2 = fig.add_subplot(2, 3, 2)
        ax2.hist(displacement*1000, bins=20, density=True, alpha=0.7, color='lightcoral', edgecolor='black')
        ax2.axvline(np.mean(displacement)*1000, color='red', linestyle='--', linewidth=2, label='均值')
        ax2.set_xlabel('位移 (mm)')
        ax2.set_ylabel('概率密度')
        ax2.set_title('输出分布')
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        
        # 3. 散点图
        ax3 = fig.add_subplot(2, 3, 3)
        ax3.scatter(E_samples/1e9, displacement*1000, alpha=0.5, s=20)
        ax3.set_xlabel('弹性模量 (GPa)')
        ax3.set_ylabel('位移 (mm)')
        ax3.set_title('输入-输出关系')
        ax3.grid(True, alpha=0.3)
        
        # 4. 收敛曲线
        ax4 = fig.add_subplot(2, 3, 4)
        cumulative_mean = np.cumsum(displacement) / np.arange(1, n_samples+1)
        ax4.plot(range(1, n_samples+1), cumulative_mean*1000, linewidth=1.5, color='blue')
        ax4.set_xlabel('样本数')
        ax4.set_ylabel('累积均值 (mm)')
        ax4.set_title('均值收敛')
        ax4.grid(True, alpha=0.3)
        
        # 5. 置信区间演化
        ax5 = fig.add_subplot(2, 3, 5)
        if n_samples >= 10:
            ci_lower = [np.percentile(displacement[:i], 2.5) for i in range(10, n_samples+1)]
            ci_upper = [np.percentile(displacement[:i], 97.5) for i in range(10, n_samples+1)]
            x_range = range(10, n_samples+1)
            ax5.fill_between(x_range, np.array(ci_lower)*1000, np.array(ci_upper)*1000, 
                            alpha=0.3, color='green', label='95% CI')
            ax5.plot(x_range, np.array(ci_lower)*1000, 'g--', linewidth=1)
            ax5.plot(x_range, np.array(ci_upper)*1000, 'g--', linewidth=1)
        ax5.set_xlabel('样本数')
        ax5.set_ylabel('位移 (mm)')
        ax5.set_title('置信区间收敛')
        ax5.legend()
        ax5.grid(True, alpha=0.3)
        
        # 6. 统计量表格
        ax6 = fig.add_subplot(2, 3, 6)
        ax6.axis('off')
        
        stats_text = f"""
        统计摘要 (n={n_samples})
        
        均值: {np.mean(displacement)*1000:.4f} mm
        标准差: {np.std(displacement)*1000:.4f} mm
        最小值: {np.min(displacement)*1000:.4f} mm
        最大值: {np.max(displacement)*1000:.4f} mm
        
        95%置信区间:
        [{np.percentile(displacement, 2.5)*1000:.4f}, 
         {np.percentile(displacement, 97.5)*1000:.4f}] mm
        """
        ax6.text(0.1, 0.5, stats_text, fontsize=12, verticalalignment='center',
                family='monospace', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
        
        plt.tight_layout()
        
        temp_file = f'unc_frame_{frame_idx:03d}.png'
        plt.savefig(temp_file, dpi=100, bbox_inches='tight')
        temp_files.append(temp_file)
        plt.close()
    
    # 生成GIF
    from PIL import Image
    images = []
    target_size = None
    
    for temp_file in temp_files:
        img = Image.open(temp_file)
        if target_size is None:
            target_size = img.size
        else:
            img = img.resize(target_size, Image.Resampling.LANCZOS)
        images.append(np.array(img))
        os.remove(temp_file)
    
    imageio.mimsave(save_path, images, fps=2)
    print(f'动画生成完成: {save_path}')

# 主程序
if __name__ == '__main__':
    print('=' * 70)
    print('结构健康监测中的不确定性量化技术')
    print('=' * 70)
    
    # 1. 蒙特卡洛模拟
    print('\n1. 蒙特卡洛模拟')
    displacement, mean_disp, std_disp, ci = monte_carlo_simulation(n_samples=10000)
    
    # 2. 多项式混沌展开
    print('\n2. 多项式混沌展开')
    coeffs, xi_train, disp_train = polynomial_chaos_expansion(order=3, n_samples=1000)
    
    # 3. 区间分析
    print('\n3. 区间分析')
    disp_interval, E_interval, F_interval = interval_analysis()
    
    # 4. 贝叶斯推断
    print('\n4. 贝叶斯推断')
    E_grid, posteriors, observations, stats = bayesian_inference(n_observations=50)
    
    # 5. 敏感性分析
    print('\n5. 敏感性分析')
    sobol_indices, total_var = sensitivity_analysis()
    
    # 6. 可靠性分析
    print('\n6. 可靠性分析')
    Pf, beta, g_values = reliability_analysis()
    
    # 7. 可视化
    print('\n7. 生成可视化...')
    visualize_uncertainty_analysis()
    
    # 8. 创建动画
    print('\n8. 生成动画...')
    create_uncertainty_animation()
    
    print('\n' + '=' * 70)
    print('仿真完成！')
    print('=' * 70)
    
    # 总结
    print('\n不确定性量化结果总结:')
    print(f'1. 蒙特卡洛: 均值={mean_disp:.6f} m, 标准差={std_disp:.6f} m')
    print(f'2. PCE: 均值={coeffs[0]:.6f} m')
    print(f'3. 区间分析: [{disp_interval[0]:.6f}, {disp_interval[1]:.6f}] m')
    print(f'4. 贝叶斯: MAP={stats[0]/1e9:.2f} GPa')
    print(f'5. 敏感性: E={sobol_indices["E"]:.3f}, F={sobol_indices["F"]:.3f}')
    print(f'6. 可靠性: Pf={Pf:.6f}, β={beta:.2f}')

---

案例研究

8.1 案例一：大跨桥梁模态参数识别的不确定性

背景
某大跨斜拉桥进行了环境振动测试，需要评估识别模态参数的不确定性。

方法

• 使用随机子空间识别（SSI）方法识别模态参数
• 采用Bootstrap重采样估计统计不确定性
• 使用贝叶斯方法量化参数后验分布
• 进行模型选择不确定性分析

结果

• 第一阶频率：0.25 ± 0.01 Hz（95%置信区间）
• 阻尼比：1.2% ± 0.3%（变异系数25%）
• 振型置信度：MAC > 0.95

关键发现

• 低阶模态识别精度高，高阶模态不确定性大
• 阻尼比识别不确定性显著高于频率
• 环境温度变化是主要不确定性来源

8.2 案例二：风力发电机塔架疲劳寿命预测

背景
某海上风力发电机需要评估塔架在随机风荷载作用下的疲劳寿命。

方法

• 建立塔架有限元模型
• 使用概率S-N曲线描述材料疲劳特性
• 采用雨流计数法提取疲劳载荷谱
• 使用蒙特卡洛模拟预测寿命分布

结果

• 设计寿命：20年
• 预测寿命中值：28年
• 5%分位寿命：18年
• 失效概率（20年内）：12%

关键发现

• 焊接细节是疲劳寿命的关键
• 风荷载不确定性对寿命预测影响最大
• 腐蚀效应显著降低疲劳寿命

8.3 案例三：核电站管道裂纹检测的不确定性

背景
某核电站管道需要评估超声检测对裂纹识别的可靠性。

方法

• 建立超声检测概率模型（POD曲线）
• 考虑检测人员、设备、环境的不确定性
• 使用贝叶斯更新融合多次检测结果
• 进行风险告知的决策分析

结果

• 检测阈值：2mm裂纹深度
• 检出概率（2mm）：70%
• 虚警概率：5%
• 综合检测可靠性：85%

关键发现

• 检测人员经验显著影响检出率
• 表面粗糙度降低检测灵敏度
• 多次检测可显著提高可靠性

---