392 工程科学学报，第42卷，第3期 ufdo<,u)为40的FFT变换 式中：Cp=JRm 较小的阀值以便更好地保护有效信号；当某一分 小波变换处理数据时基函数的选择尤为重 解方向中剪切系数值较小时，则该方向的系数主 要，小波波形越接近待处理的瞬态信号波形，处理 要是噪声系数，应该采取较大的阀值以便能压制 效果越理想.综合时频域的分辨率来看，DB族小 更多的噪声.因此，在经典剪切算法中6，仅设置 波是比较适合分析探地雷达信号的一种小波基 一个随尺度变化的阀值难以满足方向上能量变化 函数啊 的需求，容易产生“过扼杀”、“除不净”的现象 1.1.2剪切变换基本原理 本文在经典尺度阀值函数的基础上，根据有 剪切变换又称剪切波变换，最初是由Guo和 效信号和干扰信号在剪切域中不同尺度、不同方 Easley等l&-l9通过对剪切基函数进行缩放、剪切 向上能量的差异，设置一个随尺度和方向变化的 和平移合成具有膨胀性的仿射系统.当维数为 阀值函数，从而使得在剪切域进行噪声压制时每 2时，该合成膨胀仿射系统定义为： 一个子带都可以根据其能量特征自适应选择最优 阀值，其表达式如下： pAB(）={j(m)=ldetA(BAx-k):leZ,k∈Z2 (3) Thin=oV2logN×2-J-D cos 式中：山∈L2(R2),A、B为2阶的可逆矩阵且 ldet B=l,矩阵A控制尺度参数，可进行伸缩放大 j=1,2,J;1=1,2,.Li 操作：矩阵B控制着方向参数，可进行旋转剪切 (7) 操作 式中：σ为噪声标准差，N为信号长度，ε为任意小 当中AB(W)满足Parseval框架时，称该系统为复 的正数，L为尺度j上方向分解总数，e表示尺度 合小波系统.当A=A,B=B,时，复合小波系统进 j方向1上系数的能量，max(e)表示尺度j上各分 一步简化为剪切波系统.其中 解方向中系数的最大能量 A0=(各9)B=(01) 基于剪切变换的自适应阀值去噪具体过程包 括以下步骤： 对于已知信号f∈L2(R),其连续剪切变换定 (1)信号常规处理，包括：直达波去除、直流去 义为： 漂移、信号增益等； Sril,k=〈f,),(,k)∈S (4) (2)格式转换，将雷达数据格式转换为矩阵形 式中：j>0是尺度参数，s∈R是剪切参数，k∈R2是 式以保证算法的实现： 位置参数 (3)伪极坐标变换，将数据从笛卡尔坐标系 按照以下方式采样，可将剪切变换离散化： 转换到伪极坐标系，并且在伪极坐标系上进行 aj=22j,jeZ:sk=k2,-2j≤1≤2 (5) 多尺度划分，并生成剪切基函数对数据进行窗口 离散剪切系统为： 子带化： {w9(w=232w0(B'Ax-k):j≥0， (4)计算剪切系数，得到系数矩阵之后通过上 述自适应阀值函数对剪切系数矩阵进行干扰压制 -2J≤1≤2J-1,k∈Z2 (6) 处理； 12自适应阀值去噪 (5)根据去噪处理后的剪切系数对探地雷达 含噪雷达信号经过剪切变换后，数据信息被 二维数据进行重构 映射到不同尺度不同方向的剪切系数上，通常有 2 随机干扰消除 效信号会根据自身特点集中在一些特定方向上， 而随机噪声信号不具有方向性，有效信号所映射 2.1正演模拟及干扰设置 的某些特定方向上的剪切系数值往往较大，而随 利用时域有限差分法分别模拟圆形空洞与方 机噪声将在各个方向上分布，其对应的剪切系数 形空洞两种异常体，得到纯净的雷达数据.图1和 值往往较小. 图2分别为异常体几何模型和波形堆积图，从图中 从阀值去噪角度来讲，如果某一分解方向中 可以看出异常体反射界面清晰，同相轴特征明显 剪切系数值较大，则该方向为有效信号系数主要 在正演模拟获得的纯净数据中加入高斯白噪 集中的方向，对于该方向上的系数处理应该采取 声，加噪后图像信噪比为-2.5dB左右，波形堆积Cφ = r R |Ψ(ω)| 2 |ω| 式中： dω < ∞，Ψ(ω) 为 φ(t) 的 FFT 变换. 小波变换处理数据时基函数的选择尤为重 要，小波波形越接近待处理的瞬态信号波形，处理 效果越理想. 综合时频域的分辨率来看，DB 族小 波是比较适合分析探地雷达信号的一种小波基 函数[15] . 1.1.2    剪切变换基本原理 剪切变换又称剪切波变换，最初是由 Guo 和 Easley 等[18‒ 19] 通过对剪切基函数进行缩放、剪切 和平移合成具有膨胀性的仿射系统. 当维数为 2 时，该合成膨胀仿射系统定义为： ϕAB(ψ) = { ψj,l,k(x) = |detA| j 2 ψ(B lA j x−k) : j,l ∈ Z, k ∈ Z 2 } （3） ψ ∈ L 2 ( R 2 ) |detB| = 1 式 中 ： ， A、 B 为 2 阶 的 可 逆 矩 阵 且 ，矩阵 A 控制尺度参数，可进行伸缩放大 操作；矩阵 B 控制着方向参数，可进行旋转剪切 操作. 当 ϕAB (ψ) 满足 Parseval 框架时，称该系统为复 合小波系统. 当 A=A0，B=B0 时，复合小波系统进 一步简化为剪切波系统. 其中 A0 = ( 4 0 0 2 ) , B0 = ( 1 1 0 1 ) . f ∈ L 2 对于已知信号 (R) ，其连续剪切变换定 义为： S f(j,l, k) = ⟨ f,ψj,l,k ⟩ ,(j,l, k) ∈ S （4） j > 0 s ∈ R k ∈ R 式中： 是尺度参数， 是剪切参数， 2 是 位置参数. 按照以下方式采样，可将剪切变换离散化： aj = 2 −2 j , j ∈ Z;sj,k = k2 −j ,−2 −j ⩽ l ⩽ 2 −j （5） 离散剪切系统为： { ψ (0) j,l,k (x) = 2 3 j/2ψ (0) (B lA j x−k) : j ⩾ 0, −2 j ⩽ l ⩽ 2 j −1, k ∈ Z 2 } （6） 1.2    自适应阀值去噪 含噪雷达信号经过剪切变换后，数据信息被 映射到不同尺度不同方向的剪切系数上，通常有 效信号会根据自身特点集中在一些特定方向上， 而随机噪声信号不具有方向性. 有效信号所映射 的某些特定方向上的剪切系数值往往较大，而随 机噪声将在各个方向上分布，其对应的剪切系数 值往往较小. 从阀值去噪角度来讲，如果某一分解方向中 剪切系数值较大，则该方向为有效信号系数主要 集中的方向，对于该方向上的系数处理应该采取 较小的阀值以便更好地保护有效信号；当某一分 解方向中剪切系数值较小时，则该方向的系数主 要是噪声系数，应该采取较大的阀值以便能压制 更多的噪声. 因此，在经典剪切算法中[16] ，仅设置 一个随尺度变化的阀值难以满足方向上能量变化 的需求，容易产生“过扼杀”、“除不净”的现象. 本文在经典尺度阀值函数的基础上，根据有 效信号和干扰信号在剪切域中不同尺度、不同方 向上能量的差异，设置一个随尺度和方向变化的 阀值函数，从而使得在剪切域进行噪声压制时每 一个子带都可以根据其能量特征自适应选择最优 阀值，其表达式如下： Thj,l =σ √ 2logN ×2 −(J−j) ×   cos ej,l max 1⩽l⩽Lj ( ej,l ) +ε   , j = 1,2,...J;l = 1,2,...Lj （7） σ ε Lj ej,l max(e j,l) 式中： 为噪声标准差，N 为信号长度， 为任意小 的正数， 为尺度 j 上方向分解总数， 表示尺度 j 方向 l 上系数的能量， 表示尺度 j 上各分 解方向中系数的最大能量. 基于剪切变换的自适应阀值去噪具体过程包 括以下步骤： （1）信号常规处理，包括：直达波去除、直流去 漂移、信号增益等； （2）格式转换，将雷达数据格式转换为矩阵形 式以保证算法的实现； （ 3）伪极坐标变换，将数据从笛卡尔坐标系 转换到伪极坐标系，并且在伪极坐标系上进行 多尺度划分，并生成剪切基函数对数据进行窗口 子带化； （4）计算剪切系数，得到系数矩阵之后通过上 述自适应阀值函数对剪切系数矩阵进行干扰压制 处理； （5）根据去噪处理后的剪切系数对探地雷达 二维数据进行重构. 2    随机干扰消除 2.1    正演模拟及干扰设置 利用时域有限差分法分别模拟圆形空洞与方 形空洞两种异常体，得到纯净的雷达数据. 图 1 和 图 2 分别为异常体几何模型和波形堆积图，从图中 可以看出异常体反射界面清晰，同相轴特征明显. 在正演模拟获得的纯净数据中加入高斯白噪 声，加噪后图像信噪比为‒2.5 dB 左右，波形堆积 · 392 · 工程科学学报，第 42 卷，第 3 期