结构健康监测仿真-主题041-结构健康监测中的声发射技术

引言

声发射（Acoustic Emission，AE）技术是一种动态无损检测方法，通过接收材料在受力过程中释放的弹性波来监测结构内部的活动。这种技术能够实时捕捉材料内部裂纹萌生、扩展，纤维断裂，以及界面脱粘等微观损伤过程，被誉为"结构健康的听诊器"。

声发射现象最早由德国人Kaiser在1950年代系统研究，他发现金属材料在受力过程中会发出可检测的声波信号。经过70多年的发展，声发射技术已成为结构健康监测领域最重要的手段之一，广泛应用于石油化工、航空航天、电力能源、交通运输等关键基础设施的安全监测。

本章将系统介绍声发射技术的基本原理、信号处理方法、源定位技术及其在结构健康监测中的应用。

1. 声发射技术基础

1.1 声发射现象

声发射是指材料在受外力或内力作用下，由于内部结构变化（如裂纹扩展、位错运动、相变等）而释放应变能，以弹性波形式向外传播的现象。

声发射的物理机制：

1. 裂纹萌生和扩展：新裂纹面的形成和现有裂纹的扩展会释放大量能量
2. 位错运动：金属晶体中的位错滑移和增殖产生声发射信号
3. 纤维断裂：复合材料中纤维的断裂产生高频声发射信号
4. 界面脱粘：纤维与基体界面的分离产生特征声发射信号
5. 相变：材料相变过程中的体积变化产生声发射
6. 摩擦：裂纹面之间的摩擦产生连续型声发射信号

声发射信号的特点：

• 动态性：声发射是材料变形过程中的伴生现象，具有实时性
• 不可逆性：凯塞效应（Kaiser Effect）表明，材料在应力未达到历史最高水平前不会产生声发射
• 敏感性：能够检测到微米级的裂纹扩展
• 全范围监测：一个传感器可以监测数米范围内的结构活动

1.2 声发射波的类型

声发射波在固体介质中传播时，会产生多种波型：

体波：

1. 纵波（P波）：质点振动方向与波的传播方向平行，传播速度最快
   $$v_p=\sqrt{\frac{E(1-\nu )}{\rho (1+\nu )(1-2\nu )}}$$

2. 横波（S波）：质点振动方向与波的传播方向垂直，传播速度较慢
   $$v_s=\sqrt{\frac{E}{2\rho (1+\nu )}}$$

其中$E$为弹性模量，$\nu$为泊松比，$\rho$为密度。

表面波（瑞利波）：
沿材料表面传播，能量集中在表面附近一个波长范围内，是板状结构中最主要的声发射波型：
$$v_r\approx \frac{0.862+1.14\nu }{1+\nu }{v}_s$$

板波（兰姆波）：
在薄板结构中，由于上下表面的反射和干涉，形成复杂的板波模式，包括对称模式（S0, S1, S2…）和反对称模式（A0, A1, A2…）。

波型与频率的关系：

• 高频信号（>500kHz）：主要由裂纹扩展、纤维断裂产生，衰减快
• 中频信号（100-500kHz）：综合信息，传播距离适中
• 低频信号（<100kHz）：传播距离远，但分辨率低

1.3 声发射信号特征

声发射信号可以用多种参数描述：

时域参数：

1. 撞击计数（Hit Count）：超过阈值的信号个数
2. 振铃计数（Ring-down Count）：信号超过阈值的振荡次数
3. 能量（Energy）：信号幅值的平方积分
   $$E={\int }_0^T{V}^2(t)dt$$
4. 幅度（Amplitude）：信号峰值，单位dB
   $$A_{dB}=20{\log }_{10}\left(\frac{V_{peak}}{V_{ref}}\right)$$
5. 持续时间（Duration）：信号超过阈值的时间长度
6. 上升时间（Rise Time）：从信号起始到峰值的时间

频域参数：

1. 峰值频率（Peak Frequency）：功率谱密度最大的频率
2. 中心频率（Centroid Frequency）：频率的加权平均
   $$f_c=\frac{{\int }_0^{\infty }f\cdot P(f)df}{{\int }_0^{\infty }P(f)df}$$
3. 频率质心（Frequency Centroid）：与中心频率类似
4. 带宽（Bandwidth）：信号频率分布范围

统计参数：

1. b值：描述声发射事件幅度分布的参数，与裂纹扩展机制相关
2. RA值（Rise Time/Amplitude）：上升时间与幅度的比值
3. 平均频率（Average Frequency）：振铃计数与持续时间的比值

1.4 凯塞效应与费利西蒂比

凯塞效应（Kaiser Effect）：
材料在单调加载过程中，当应力未超过历史最大应力时，声发射活动微弱或不产生声发射。这一现象由Kaiser于1950年发现，是声发射技术的重要理论基础。

凯塞效应的物理本质：

• 材料内部的微缺陷在首次加载时已经激活
• 再次加载时，只有应力超过历史最高水平，新的损伤才会产生
• 适用于大多数金属材料和复合材料

费利西蒂比（Felicity Ratio）：
定义为产生明显声发射的应力与历史最大应力的比值：
$$FR=\frac{{\sigma }_{AE}}{{\sigma }_{max}}$$

• $FR\ge 0.95$：结构状态良好
• $0.80\le FR<0.95$：需要关注
• $FR<0.80$：存在严重损伤

费利西蒂比突破了凯塞效应的限制，可用于评估结构的损伤程度。

2. 声发射检测系统

2.1 系统组成

声发射检测系统主要由以下部分组成：

传感器（换能器）：

• 压电式传感器：最常用，利用压电效应将机械振动转换为电信号
• 电容式传感器：灵敏度高，但易受环境影响
• 光纤传感器：抗电磁干扰，适合恶劣环境
• MEMS传感器：微型化，适合密集阵列

传感器的关键参数：

• 谐振频率：通常在150kHz-400kHz之间
• 灵敏度：单位声压产生的电压输出
• 频率响应：有效工作频率范围

前置放大器：

• 增益：通常40dB或60dB
• 带宽：与传感器匹配
• 噪声系数：影响系统最小检测幅度

信号采集系统：

• 阈值设定：过滤背景噪声
• 采样率：通常1-10MHz
• 分辨率：12-16位ADC
• 通道数：从单通道到数百通道

分析软件：

• 实时参数提取
• 源定位计算
• 模式识别
• 数据存储与管理

2.2 传感器布置

布置原则：

1. 覆盖范围：根据材料衰减特性确定传感器间距

   • 金属：3-5米
   • 混凝土：1-2米
   • 复合材料：0.5-1米

2. 几何配置：

   • 线性布置：适用于管道、钢索
   • 平面布置：适用于板壳结构
   • 立体布置：适用于复杂结构

3. 定位精度：

   • 传感器间距越小，定位精度越高
   • 至少需要3个传感器进行平面定位
   • 至少需要4个传感器进行空间定位

常用布置方式：

三角形布置：
三个传感器构成等边三角形，适合圆形或三角形区域监测。

正方形布置：
四个传感器构成正方形，适合矩形区域监测，定位精度高。

阵列布置：
多个传感器按规则网格布置，适合大面积监测。

2.3 噪声控制

噪声来源：

1. 机械噪声：摩擦、碰撞、流体流动
2. 电磁噪声：电源干扰、无线电干扰
3. 环境噪声：风雨、温度变化
4. 系统噪声：电子元件热噪声

噪声抑制方法：

硬件方法：

• 屏蔽电缆和接地
• 带通滤波器
• 差分放大
• 屏蔽罩

软件方法：

• 阈值设定
• 时间门控
• 模式识别
• 小波去噪

现场措施：

• 减少机械接触
• 控制环境条件
• 选择合适的工作频段

3. 声发射信号处理

3.1 时域分析

参数分析法：
提取信号的基本特征参数，如幅度、能量、计数、持续时间等。这种方法计算简单，适合实时监测。

波形分析法：
直接分析信号的时域波形，可以获取更详细的信息：

• 波形形态特征
• 到达时间精确测量
• 波形相关分析

统计分析方法：

1. 累积计数分析：
   观察声发射活动随时间或载荷的变化趋势

2. 幅度分布分析：
   分析不同幅度事件的分布规律
   $$N(A)=N_0\cdot A^{-b}$$
   其中$b$值反映损伤机制：

   • $b>1.5$：以微裂纹为主
   • $b<1.5$：以宏观裂纹为主

3. 事件间隔分析：
   分析连续事件之间的时间间隔分布

3.2 频域分析

快速傅里叶变换（FFT）：
将时域信号转换到频域，分析频率成分：
$$X(f)={\int }_{-\infty }^{\infty }x(t)e^{-j2\pi ft}dt$$

功率谱密度（PSD）：
$$P(f)=|X(f){|}^2$$

频域特征：

• 峰值频率：主要频率成分
• 频带能量：不同频段的能量分布
• 频谱质心：频率分布的中心位置

频域分析的应用：

• 识别不同损伤机制（纤维断裂vs基体开裂）
• 滤波器设计
• 传播特性研究

3.3 时频分析

短时傅里叶变换（STFT）：
$$STFT(t,f)={\int }_{-\infty }^{\infty }x(\tau )w(\tau -t)e^{-j2\pi f\tau }d\tau$$

优点：直观显示频率随时间的变化
缺点：时间分辨率和频率分辨率相互制约

小波变换（Wavelet Transform）：
$$WT(a,b)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{\infty }x(t){\psi }^{\ast }\left(\frac{t-b}{a}\right)dt$$

优点：多分辨率分析，适合非平稳信号
常用小波：Morlet小波、Daubechies小波

希尔伯特-黄变换（HHT）：

1. 经验模态分解（EMD）将信号分解为IMF分量
2. 对每个IMF进行希尔伯特变换得到瞬时频率

优点：自适应分解，适合非线性非平稳信号

3.4 模式识别

特征提取：
从声发射信号中提取能够区分不同损伤模式的特征：

• 时域特征：幅度、能量、持续时间、RA值
• 频域特征：峰值频率、中心频率、带宽
• 时频特征：小波能量分布

分类方法：

1. 无监督学习：

   • K-means聚类
   • 层次聚类
   • 自组织映射（SOM）

2. 有监督学习：

   • 支持向量机（SVM）
   • 人工神经网络（ANN）
   • 随机森林
   • 深度学习

损伤模式分类：

• 基体开裂
• 纤维断裂
• 界面脱粘
• 分层
• 摩擦

4. 声发射源定位

4.1 时差定位法

基本原理：
利用声发射信号到达不同传感器的时间差，通过求解双曲线方程组确定源位置。

平面定位（二维）：
假设波速$v$已知，对于传感器$i$和$j$：
$$\sqrt{(x-x_i{)}^2+(y-y_i{)}^2}-\sqrt{(x-x_j{)}^2+(y-y_j{)}^2}=v\cdot \Delta t_{ij}$$

其中$(x,y)$为源坐标，$(x_i,y_i)$为传感器坐标，$\Delta t_{ij}$为到达时间差。

求解方法：

1. 解析法：对于3个传感器，可以直接求解
2. 迭代法：对于更多传感器，使用最小二乘法优化
3. 网格搜索法：在区域内搜索使误差最小的点

定位误差来源：

• 波速不确定性
• 到达时间测量误差
• 传感器位置误差
• 各向异性材料

4.2 幅值衰减定位法

基本原理：
利用声发射信号幅度随传播距离衰减的特性进行定位。

衰减模型：
$$A(r)=A_0\cdot r^{-\alpha }\cdot e^{-\beta r}$$

其中：

• $A_0$：源处幅度
• $r$：传播距离
• $\alpha$：几何衰减指数（通常为1）
• $\beta$：材料衰减系数

定位方法：
对于多个传感器，建立方程组：
$$\frac{A_i}{A_j}={\left(\frac{r_j}{r_i}\right)}^{\alpha }\cdot e^{-\beta (r_i-r_j)}$$

通过优化求解源位置。

适用条件：

• 材料衰减特性已知
• 幅度测量准确
• 传播路径简单

4.3 三角测量法

基本原理：
利用多个传感器接收到的信号，通过几何关系确定源位置。

三维定位：
至少需要4个不共面的传感器，求解以下方程组：
$$\sqrt{(x-x_i{)}^2+(y-y_i{)}^2+(z-z_i{)}^2}=v\cdot t_i$$

其中$t_i$为信号到达传感器$i$的时间。

优化目标函数：
$$\underset{x,y,z}{\min }\sum _{i=1}^n{\left(\sqrt{(x-x_i{)}^2+(y-y_i{)}^2+(z-z_i{)}^2}-v\cdot t_i\right)}^2$$

4.4 先进定位技术

基于波形的定位：
利用完整波形信息进行定位，考虑波的传播特性：

• 波型识别（纵波、横波、表面波）
• 模态分析
• 频散补偿

概率定位法：
考虑各种不确定性因素，给出源位置的概率分布：
$$P(x,y,z)=\prod _{i=1}^n{P}_i(x,y,z)$$

机器学习定位：
训练神经网络直接预测源位置：

• 输入：各传感器的特征参数
• 输出：源坐标
• 优点：不需要精确的波速信息

5. 声发射在结构健康监测中的应用

5.1 金属材料监测

裂纹监测：

• 检测疲劳裂纹的萌生和扩展
• 评估裂纹扩展速率
• 预测剩余寿命

压力容器监测：

• 检测腐蚀、应力腐蚀开裂
• 评估焊接缺陷
• 水压试验监测

管道监测：

• 检测腐蚀、裂纹、泄漏
• 评估管道完整性
• 定位损伤位置

典型案例：
某石化企业加氢反应器声发射监测：

• 在役监测期间检测到多处声发射源
• 定位分析发现主要集中在焊缝区域
• 停机检查确认了裂纹的存在
• 避免了可能的灾难性事故

5.2 复合材料监测

损伤机制识别：
通过声发射参数识别不同损伤类型：

• 基体开裂：低幅度、高频率
• 纤维断裂：高幅度、低频率
• 界面脱粘：中等幅度、中等频率
• 分层：持续时间长、频率低

损伤演化监测：

• 记录损伤累积过程
• 评估损伤严重程度
• 预测剩余强度

典型案例：
风力发电机叶片声发射监测：

• 在疲劳试验中实时监测损伤演化
• 识别出基体开裂→界面脱粘→纤维断裂的损伤序列
• 建立了损伤与声发射参数的定量关系
• 为叶片设计优化提供了依据

5.3 混凝土结构监测

裂缝监测：

• 检测新裂缝的产生
• 监测现有裂缝的扩展
• 评估裂缝活动性

预应力损失监测：

• 检测预应力筋的腐蚀和断裂
• 评估预应力损失程度
• 预测结构承载力

典型案例：
桥梁混凝土箱梁声发射监测：

• 长期监测发现裂缝活动规律
• 识别出车辆荷载引起的周期性声发射
• 评估了裂缝的稳定性
• 为维修决策提供了依据

5.4 岩土工程监测

边坡稳定性监测：

• 检测岩体微破裂
• 评估边坡稳定性
• 预警滑坡灾害

隧道监测：

• 检测围岩松动
• 监测支护结构损伤
• 评估施工安全

矿山监测：

• 检测岩爆前兆
• 监测采空区稳定性
• 预警地质灾害

5.5 其他应用

轴承监测：

• 检测疲劳剥落
• 评估润滑状态
• 预测剩余寿命

刀具磨损监测：

• 检测刀具破损
• 评估磨损程度
• 优化加工参数

泄漏检测：

• 检测管道、容器泄漏
• 定位泄漏点
• 评估泄漏速率

6. Python声发射仿真实现

6.1 声发射信号生成

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft, fftfreq

def generate_ae_signal(duration=0.001, sample_rate=1000000, 
                       frequency=150000, amplitude=1.0, 
                       decay_time=0.0002, noise_level=0.05):
    """
    生成合成声发射信号
    
    Parameters:
    -----------
    duration : float
        信号持续时间（秒）
    sample_rate : int
        采样率（Hz）
    frequency : float
        中心频率（Hz）
    amplitude : float
        信号幅度
    decay_time : float
        衰减时间常数（秒）
    noise_level : float
        噪声水平
    
    Returns:
    --------
    t : ndarray
        时间向量
    signal : ndarray
        声发射信号
    """
    # 时间向量
    t = np.linspace(0, duration, int(sample_rate * duration))
    
    # 包络（指数衰减）
    envelope = np.exp(-t / decay_time)
    
    # 载波信号（调制）
    carrier = np.sin(2 * np.pi * frequency * t)
    
    # 调制信号
    ae_signal = amplitude * envelope * carrier
    
    # 添加高斯噪声
    noise = noise_level * amplitude * np.random.randn(len(t))
    ae_signal += noise
    
    return t, ae_signal

def generate_burst_signal(t, arrival_time, frequency=150000, 
                          amplitude=1.0, decay_time=0.0001):
    """
    生成突发型声发射信号
    
    Parameters:
    -----------
    t : ndarray
        时间向量
    arrival_time : float
        到达时间（秒）
    frequency : float
        频率（Hz）
    amplitude : float
        幅度
    decay_time : float
        衰减时间（秒）
    
    Returns:
    --------
    signal : ndarray
        突发信号
    """
    signal = np.zeros_like(t)
    
    # 找到到达时间对应的索引
    idx = np.argmin(np.abs(t - arrival_time))
    
    if idx < len(t):
        # 生成衰减信号
        t_local = t[idx:] - t[idx]
        envelope = np.exp(-t_local / decay_time)
        carrier = np.sin(2 * np.pi * frequency * t_local)
        burst = amplitude * envelope * carrier
        
        # 填充到信号中
        signal[idx:idx+len(burst)] = burst[:len(t)-idx]
    
    return signal

# 生成示例信号
if __name__ == "__main__":
    # 参数设置
    sample_rate = 2000000  # 2MHz
    duration = 0.005  # 5ms
    
    # 生成突发型声发射信号
    t = np.linspace(0, duration, int(sample_rate * duration))
    
    # 模拟3个声发射事件
    ae1 = generate_burst_signal(t, 0.001, frequency=150000, amplitude=1.0)
    ae2 = generate_burst_signal(t, 0.0025, frequency=300000, amplitude=0.6)
    ae3 = generate_burst_signal(t, 0.0038, frequency=200000, amplitude=0.8)
    
    # 合成信号
    ae_signal = ae1 + ae2 + ae3 + 0.02 * np.random.randn(len(t))
    
    # 绘图
    fig, axes = plt.subplots(3, 1, figsize=(12, 8))
    
    # 时域波形
    axes[0].plot(t * 1000, ae_signal)
    axes[0].set_xlabel('Time (ms)')
    axes[0].set_ylabel('Amplitude')
    axes[0].set_title('Synthetic AE Signal')
    axes[0].grid(True, alpha=0.3)
    
    # 频谱
    freqs = fftfreq(len(t), 1/sample_rate)
    spectrum = np.abs(fft(ae_signal))
    
    axes[1].plot(freqs[:len(freqs)//2] / 1000, spectrum[:len(spectrum)//2])
    axes[1].set_xlabel('Frequency (kHz)')
    axes[1].set_ylabel('Magnitude')
    axes[1].set_title('Frequency Spectrum')
    axes[1].set_xlim(0, 500)
    axes[1].grid(True, alpha=0.3)
    
    # 短时傅里叶变换
    f, t_stft, Zxx = signal.stft(ae_signal, sample_rate, nperseg=1024)
    
    im = axes[2].pcolormesh(t_stft * 1000, f / 1000, np.abs(Zxx), 
                            shading='gouraud', cmap='jet')
    axes[2].set_xlabel('Time (ms)')
    axes[2].set_ylabel('Frequency (kHz)')
    axes[2].set_title('STFT Spectrogram')
    axes[2].set_ylim(0, 500)
    plt.colorbar(im, ax=axes[2], label='Magnitude')
    
    plt.tight_layout()
    plt.savefig('ae_signal_example.png', dpi=150)
    plt.show()

6.2 声发射参数提取

import numpy as np
from scipy.signal import find_peaks

class AEParameterExtractor:
    """声发射参数提取器"""
    
    def __init__(self, sample_rate=1000000, threshold=0.1):
        """
        初始化参数提取器
        
        Parameters:
        -----------
        sample_rate : int
            采样率（Hz）
        threshold : float
            检测阈值
        """
        self.sample_rate = sample_rate
        self.threshold = threshold
        self.dt = 1.0 / sample_rate
        
    def extract_hit(self, signal, start_idx, end_idx):
        """
        提取单个撞击的参数
        
        Parameters:
        -----------
        signal : ndarray
            信号
        start_idx : int
            起始索引
        end_idx : int
            结束索引
        
        Returns:
        --------
        params : dict
            声发射参数
        """
        hit_signal = signal[start_idx:end_idx]
        t_hit = np.arange(len(hit_signal)) * self.dt
        
        # 幅度（dB）
        peak_idx = np.argmax(np.abs(hit_signal))
        amplitude_linear = np.abs(hit_signal[peak_idx])
        amplitude_db = 20 * np.log10(amplitude_linear / 1e-6)  # 参考1μV
        
        # 上升时间
        rise_time = peak_idx * self.dt
        
        # 持续时间
        duration = (end_idx - start_idx) * self.dt
        
        # 能量（平方积分）
        energy = np.sum(hit_signal**2) * self.dt
        
        # 振铃计数（过阈值次数）
        crossings = np.where(np.diff(np.signbit(hit_signal - self.threshold).astype(int)) != 0)[0]
        ring_count = len(crossings) // 2
        
        # 平均频率
        avg_frequency = ring_count / duration if duration > 0 else 0
        
        # RA值（上升时间/幅度）
        ra_value = rise_time / amplitude_linear if amplitude_linear > 0 else 0
        
        params = {
            'amplitude_db': amplitude_db,
            'amplitude_linear': amplitude_linear,
            'rise_time': rise_time * 1e6,  # 转换为微秒
            'duration': duration * 1e6,     # 转换为微秒
            'energy': energy,
            'ring_count': ring_count,
            'avg_frequency': avg_frequency / 1000,  # 转换为kHz
            'ra_value': ra_value * 1e3,  # 转换为ms/V
            'threshold': self.threshold
        }
        
        return params
    
    def process_signal(self, signal):
        """
        处理整个信号，提取所有撞击
        
        Parameters:
        -----------
        signal : ndarray
            输入信号
        
        Returns:
        --------
        hits : list
            所有撞击的参数列表
        """
        hits = []
        
        # 检测超过阈值的区域
        above_threshold = np.abs(signal) > self.threshold
        
        # 找到连续区域
        changes = np.diff(above_threshold.astype(int))
        start_indices = np.where(changes == 1)[0] + 1
        end_indices = np.where(changes == -1)[0] + 1
        
        # 处理边界情况
        if above_threshold[0]:
            start_indices = np.insert(start_indices, 0, 0)
        if above_threshold[-1]:
            end_indices = np.append(end_indices, len(signal))
        
        # 提取每个撞击的参数
        for start, end in zip(start_indices, end_indices):
            # 添加前后扩展时间
            start_ext = max(0, start - int(0.0001 * self.sample_rate))
            end_ext = min(len(signal), end + int(0.0005 * self.sample_rate))
            
            params = self.extract_hit(signal, start_ext, end_ext)
            params['arrival_time'] = start * self.dt * 1000  # 转换为ms
            hits.append(params)
        
        return hits

# 使用示例
if __name__ == "__main__":
    # 生成测试信号
    sample_rate = 1000000
    duration = 0.01
    t = np.linspace(0, duration, int(sample_rate * duration))
    
    # 生成多个声发射事件
    signal = np.zeros_like(t)
    for i in range(5):
        arrival = np.random.uniform(0.001, 0.009)
        amp = np.random.uniform(0.5, 1.0)
        freq = np.random.uniform(100000, 300000)
        signal += generate_burst_signal(t, arrival, frequency=freq, amplitude=amp)
    
    signal += 0.02 * np.random.randn(len(t))
    
    # 参数提取
    extractor = AEParameterExtractor(sample_rate=sample_rate, threshold=0.1)
    hits = extractor.process_signal(signal)
    
    print(f"检测到 {len(hits)} 个声发射事件")
    print("\n参数统计:")
    print(f"平均幅度: {np.mean([h['amplitude_db'] for h in hits]):.2f} dB")
    print(f"平均能量: {np.mean([h['energy'] for h in hits]):.6f}")
    print(f"平均持续时间: {np.mean([h['duration'] for h in hits]):.2f} μs")

6.3 声发射源定位

import numpy as np
from scipy.optimize import minimize

class AESourceLocator:
    """声发射源定位器"""
    
    def __init__(self, sensor_positions, wave_velocity):
        """
        初始化定位器
        
        Parameters:
        -----------
        sensor_positions : ndarray
            传感器位置，shape (n_sensors, 2) 或 (n_sensors, 3)
        wave_velocity : float
            波速（m/s）
        """
        self.sensor_positions = np.array(sensor_positions)
        self.wave_velocity = wave_velocity
        self.n_sensors = len(sensor_positions)
        self.dim = sensor_positions.shape[1]
        
    def locate_tdoa(self, arrival_times):
        """
        基于到达时间差（TDOA）定位
        
        Parameters:
        -----------
        arrival_times : ndarray
            各传感器的到达时间（秒）
        
        Returns:
        --------
        source_position : ndarray
            源位置
        residual : float
            残差
        """
        # 以第一个传感器为参考
        t_ref = arrival_times[0]
        time_diffs = arrival_times - t_ref
        
        # 初始猜测：传感器几何中心
        x0 = np.mean(self.sensor_positions, axis=0)
        
        # 定义目标函数
        def objective(x):
            distances = np.linalg.norm(self.sensor_positions - x, axis=1)
            predicted_diffs = (distances - distances[0]) / self.wave_velocity
            return np.sum((predicted_diffs - time_diffs)**2)
        
        # 优化
        result = minimize(objective, x0, method='Nelder-Mead')
        
        return result.x, result.fun
    
    def locate_grid_search(self, arrival_times, bounds, resolution=100):
        """
        网格搜索定位
        
        Parameters:
        -----------
        arrival_times : ndarray
            到达时间
        bounds : list
            搜索范围，如 [(xmin, xmax), (ymin, ymax)]
        resolution : int
            网格分辨率
        
        Returns:
        --------
        source_position : ndarray
            源位置
        error_map : ndarray
            误差分布
        """
        # 创建网格
        if self.dim == 2:
            x_grid = np.linspace(bounds[0][0], bounds[0][1], resolution)
            y_grid = np.linspace(bounds[1][0], bounds[1][1], resolution)
            X, Y = np.meshgrid(x_grid, y_grid)
            error_map = np.zeros_like(X)
            
            # 计算每个网格点的误差
            for i in range(resolution):
                for j in range(resolution):
                    pos = np.array([X[i, j], Y[i, j]])
                    distances = np.linalg.norm(self.sensor_positions - pos, axis=1)
                    predicted_times = distances / self.wave_velocity
                    error_map[i, j] = np.sum((predicted_times - arrival_times)**2)
            
            # 找到最小误差位置
            min_idx = np.unravel_index(np.argmin(error_map), error_map.shape)
            source_position = np.array([X[min_idx], Y[min_idx]])
            
        else:  # 3D
            # 简化处理，使用TDOA方法
            source_position, _ = self.locate_tdoa(arrival_times)
            error_map = None
        
        return source_position, error_map
    
    def calculate_location_error(self, true_position, estimated_position):
        """计算定位误差"""
        return np.linalg.norm(true_position - estimated_position)

# 定位示例
if __name__ == "__main__":
    # 传感器布置（正方形，边长1m）
    sensor_positions = np.array([
        [0, 0],
        [1, 0],
        [1, 1],
        [0, 1]
    ])
    
    # 波速（钢中，m/s）
    wave_velocity = 5000
    
    # 真实源位置
    true_source = np.array([0.6, 0.7])
    
    # 计算理论到达时间
    distances = np.linalg.norm(sensor_positions - true_source, axis=1)
    arrival_times = distances / wave_velocity
    
    # 添加噪声
    arrival_times += np.random.normal(0, 1e-6, len(arrival_times))
    
    # 定位
    locator = AESourceLocator(sensor_positions, wave_velocity)
    estimated_source, residual = locator.locate_tdoa(arrival_times)
    
    # 计算误差
    error = locator.calculate_location_error(true_source, estimated_source)
    
    print(f"真实源位置: ({true_source[0]:.3f}, {true_source[1]:.3f}) m")
    print(f"估计源位置: ({estimated_source[0]:.3f}, {estimated_source[1]:.3f}) m")
    print(f"定位误差: {error*1000:.2f} mm")
    print(f"残差: {residual:.10f}")