曲波变换在航空伽马能谱数据处理中的应用研究

李兵海，卢亚运，张光雅，李江坤，叶发旺，童勤龙，张翔

（1.核工业航测遥感中心，河北 石家庄 050002；2.核工业北京地质研究院，北京 100029）

航空伽马能谱测量是主要的铀矿勘查方法之一，除直接寻找铀矿外，还可直接寻找钍矿、钾盐矿等放射性矿产，也可应用于寻找与放射性矿产伴生的非放射性金属和非金属矿产，如金、铜、稀土矿、磷矿、油气等。航空伽马能谱测量具有不受地面地形地貌和交通条件的限制、速度快、效率高、成本低、覆盖面广、信息量大、找矿效果明显、对环境影响极小等优势，在国内外矿产勘查、远景评价及各项地质研究中占有重要的位置，起着先行作用［1］。然而，随着找矿目标层从浅到深的变化，航放数据固有的噪声严重影响了航空伽马能谱数据的应用效果，导致深部反演和解释不准确。因此，抑制复杂的噪声并提高成果数据的质量，已成为现今重要且必要的任务。

目前航空伽马能谱数据降噪主要为一维尺度，故最终结果会出现沿测线方向的条带，这种条带的处理方法一般采用微调平技术，其缺点在于一次只能调平一种宽度的条带，对于有多种宽度的复杂条带不仅需要调平多次，还容易丢失异常信息，造成了无条带和局部异常不可兼得的局面。曲波变换的出现，可以完美的解决这种情况，利用曲波多方向的特性可以降低航空伽马能谱数据中的沿测线方向的条带，利用曲波多尺度的特性可以降低航空伽马能谱数据中的白噪声，实现无条带和局部异常兼得的效果。

国内外许多学者将曲波变换运用于图像数据处理和地震资料数据处理，推动了所属行业技术的发展［2-12］。国内学者将曲波变换运用于位场数据处理方面的研究较晚，陈召曦［13］通过建立简单、复杂组合模型对曲波变换的多尺度、多方向性进行了分析， 对中国大陆及邻区卫星布格重力异常数据进行了曲波变换及阈值重构，得到了不同尺度上的重力分解结果，对中国大陆及邻区的区域重力场特征、构造分区等方面进行了分析，给出了初步的地质解释。张杨［14］等从曲波变换的基本原理入手，通过重力位场理论模型数据分析了曲波变换的多尺度分解重构能力，并且利用加噪理论模型数据分析了曲波变换的阈值降噪能力；此外，还使用曲波变换对南岭东部地区布格重力异常资料进行了有效信号提取，验证了曲波变换可同时适用于位场数据的分解和降噪处理研究，为重磁数据多尺度分析处理提供参考。

王宁、高玲琦［15-16］等将曲波变换应用于航空电磁数据处理，将曲波变换降噪与传统的降噪方法（奇异值分解、主成分分析和小波变换）做比较，得出曲波变换降噪效果最佳，重构信号最准确，信噪比最高的结论。因降噪方法与数据类型无关，故本文不再讨论与证明曲波变换降噪与传统的降噪方法的优劣。

1 基于曲波变换的降噪方法

曲波变换是在小波变换基础上发展而来，是小波变换的高维泛化，可表示不同尺度和不同方向的图像，可以提取不同尺度任意方向及任意方向组合的信息。曲波变换不仅继承了小波变换的多尺度性质，还具有多方向性，因此曲波变换能更好地表征二维或更高维的信号。开发基于曲波变换的数据条带去除和降噪方法，对有效提高航空伽马能谱数据质量，降低解释难度和成本能够起到积极的作用。基于曲波变换的调平方法只对数据的某些细节尺度和角度进行局部改变，故对数据的原始相位关系影响很小，在去除调平误差时不容易引起数据畸变。

离散曲波变换在波数域内具有紧凑支撑特性，可以将信号转换到不同的尺度和角度进行存放。我们首先在波数域定义曲波。为此，我们定义二维波数域内的径向窗W͂j(ω)和角度窗Vj(ω)为

式中：ω= (ω1，ω2)，为波数域内的变量；j 为尺度参数；［j/2］表示j 的整数部分；φ由两个一维低通窗的乘积定义，可表示为

式中：滤波器φ满足0 ≤φ≤1。由此，我们可以进一步定义如下“笛卡尔”窗，即

根据“笛卡尔”窗的定义可知，径向窗和角度窗共同构造出多个楔形支撑区间，波数域在楔形的边缘处{(ω1，ω2)|2j≤ω1≤2j+1，-2-j/2≤ω2/ω1≤2-j/2}被分割。接下来引入相同间隔的斜率tanθl=l· 2-[j/2]，其中，角度参数l满足l=-2-[j/2]，...，2[j/2]- 1，则公式（4）可写为

式中：剪切矩阵Sθ为

图1 为曲波分解5个尺度第2 尺度为16个方向时的尺度方向示意图，图中数字代表尺度，蓝线代表方向，中间标号为1 的矩形区域代表第一尺度（也称为粗尺度，低频），不具有方向特性；从第2 尺度开始，具有方向性。随着同心方环向外扩展，尺度参数增大，角度窗的数量增大，对应的信息也趋向于细节，最高尺度存储的信息频率最高，角度窗对各尺度的信息从不同方向进行划分，对信息进一步细化。一般来说中间几个尺度称为细节尺度，较高的尺度称为精细尺度。

图1 曲波分解后的尺度与方向示意图Fig. 1 Figure of scale and direction after curvelet decomposition

根据波数域中曲波的定义，我们进一步讨论空间域中曲波的形态。因此，首先引入R2内空间变量为x 的连续曲波。我们定义一个“母曲波”为φj(x)，其傅里叶变换为φ̂j(x) =Uj(ω)，其中Uj(ω)为连续曲波变换中的窗函数。引入相同间隔的旋转角度序列θl= 2πl· 2-[j/2]，l=0，1，...，以及平移参数序列k= (k1，k2) ∈Z2。通过母曲波φj(x)的旋转和平移，可以得到该尺度下所有的曲波。由此，空间域内尺度为2-j，方向为θl= 2πl· 2-[j/2]，位置为x(j，l)k= R-1θl(k1·2-j，k2· 2-j/2)的曲波可定义为

式中：旋转矩阵Rθ为

而曲波系数由φj，l，k(x) 与函数f∈L2(R2)的内积定义，可表示为

根据Plancherel 定理，曲波系数也可以在波数域内由下式得到，即

波数域与空间域内的曲波满足傅里叶变换关系，且由于傅里叶变换的90°旋转特性，使得波数域和空间域内相对应的曲波方向相互正交。

基于上述关于曲波变换理论的介绍，可以将任意航空伽马能谱数据通过曲波变换转换至一系列不同尺度、不同角度下的曲波系数。条带误差主要分布在曲波域不同尺度特定角度的系数中，对这些方向系数进行阈值处理即可有效消除条带。地质信号在曲波域内通常对应数值较大的系数，而噪声多对应数值较小的系数。因此，对曲波域各尺度不同方向系数进行阈值处理，滤除与噪声、条带相关的系数，再将剩余的曲波系数重构回空间域，即可实现航空伽马能谱数据的条带去除和降噪。

2 航空伽马能谱数据曲波变换处理方法

2.1 航空伽马能谱数据曲波分解参数设置及特征

由于曲波分解用于航空伽马能谱数据处理中的主要任务之一为消除条带，最宽的条带具有架次性，即飞行一个架次的宽度，曲波分解的尺度设置取决于条带最宽的宽度，像素最差几何分辨率有方向系数的尺度为第二尺度，第一尺度不分方向，没有方向系数，第二尺度开始具有方向系数，对于网格间距为100 m 的网格数据来说，分解为7个尺度时，每个像素代表约6.4 km 的范围，分解为6个尺度时，第二尺度每个像素代表3.2 km 的范围，目前航空伽马能谱测量的比例尺一般均为1∶5 万，线距为500 m，3.2 km 的范围一般相当于一架次测线宽度的一半， 6.4 km 大致代表12 条测线的宽度，大致可以覆盖一架次的测线，因此分解为7个尺度比较合理。

表1 曲波变换各尺度系数结构Table 1 The structure of the scale coefficients of the curvelet

分解的角度方向个数第二尺度设置为16个方向，每一方向代表22.5 度的范围，在去除条带时，范围太大，容易影响非条带数据，因此第二尺度设置为32个方向，每一方向代表11.25 度的范围，较为精准。

表2 曲波变换各尺度第一组方向矩阵形式Table 2 The matrix form of the first group of directions in each scale of the curvelet

第二尺度设置为32个方向时，曲波分解后第3、4 尺度为64个方向，第5、6 尺度为128个方向，第7 尺度为256个方向。当分解的二维数据为正方形时，曲波分解后每一个尺度的方向分为8 组，每组4个矩阵；分解的数据为长方形时，曲波分解后每一个尺度的方向分为4 组，每组8个矩阵；组与组之间矩阵形式相同，每一个方向均为一个矩阵，一组内矩阵的大小部分相同。阈值重构就是根据分解后每个方向代表的含义，选择需要处理的方向进行阈值处理，航空伽马能谱数据曲波分解后除第一尺度外，共672个方向，由于第一尺度为框架数据（也称为粗尺度，低频），不具有方向特性，无需处理，所以航空伽马能谱数据曲波变换处理就是在672个方向中选择需要处理的方向进行不同的阈值处理。条带阈值处理方法一般采用硬阈值函数，噪声阈值处理方法一般采用软阈值函数或线性阈值函数。

以某测区钍元素含量数据为例，其网格化数据为2 095×2 558 的二维数据，面积约200 km×250 km。经过上述设置的曲波变换后，可以得到系列曲波系数，其结构见表1。表2 为第二尺度至第7 尺度第一组方向矩阵形式。因为各组矩阵形式相同，所以每一尺度仅列了一组方向的矩阵形式。从矩阵形式上可以看出，随着尺度序数的增大，对应的信息也趋向于细节，对应的细节也更丰富。

2.2 不规则测区的处理

航空伽马能谱测区一般均为不规则形状，在应用其他方法（加权平均法、低通滤波法、小波法）降噪时，都会产生一定的边界效应，然而，基于曲波变换的降噪方法由于具有多尺度性，可以有效地处理不规则测区获取的数据。高玲琦［15］等用一个常数来填充缺失的区域，然后进行调平，得到了较好的效果；本文作者尝试了这个方法，发现常数数值的大小非常关键，如果给出的值过大或过小均会产生边界效应，即在边界部位产生假值（图2，填充常数为0），经过多次实验发现，当这个数值设为测区所有测量值的均值时，几乎没有边界效应。

图2 某测区钍元素含量曲波变换边界效应Fig. 2 The boundary effect map reconstructed by the curve wave transform of thorium content in the study area

2.3 实测数据处理效果

图3 为某测区利用全谱计算的钍元素含量原始数据。由图3 可见，原始数据中不仅存在随机噪声，还存在沿测线方向不同尺度（315 度方向）的条带。利用本文所述方法对该数据进行阈值处理，利用不同尺度的组合系数重构了网格数据文件，为了评价降噪效果，利用参考文献［15］中提供的均方根差（Root Mean Square Error，RMSE）和平滑度（R）公式计算了原始数据与处理后不同尺度组合重构数据的均方根差和平滑度（表3）。降噪后重构信号文件名的后两位代表了重构的尺度组合数字，例如：13为1 至3 尺度组合的阈值重构图，15 为1 至5 尺度的阈值重构图，以此类推。RMSE 体现了降噪后重构信号与原始信号之间的偏差，RMSE越小，两者偏差越小；R 体现了降噪后重构信号的平滑程度，R 值越小，数据越平滑。数据表明随着所用尺度的增加均方根差在减小，平滑度在变差。为了综合评价重构后的数据质量及性能，设计了参数加权平滑度。根据加权平滑度的数值重构后的数据质量及性能，不难发现重构数据el2019_nathw17（第1 尺度至第7 尺度的组合）平滑度R 与加权平滑度R'均有一个突变，R 增加了很多，而RMSE 没有减小很多，平滑度R 的突变预示着噪声信息的陡增。因此不建议使用1 至7 尺度组合的阈值重构数据，推荐曲波变换1 至6 尺度组合的阈值重构图作为最终提交的标准图件。

表3 不同尺度重构数据均方根误差与平滑度Table 3 The root mean square error and smoothness of reconstructed data at different scales

图3 某测区钍元素含量原始数据Fig. 3 The original data map of thorium content in study area

将曲波变换不同尺度组合阈值重构文件的加权平滑度（1 至7 尺度组合因噪声太大，不参与）绘制成曲线图（图4）。由图4 可见，1至6 尺度组合信息量最大，可为精细信息图件，1 至4 尺度组合与1 至5 尺度组合可分为一类，为细节信息图件，剩余为一类，为粗糙信息图件。

图4 不同尺度重构数据加权平滑度曲线图Fig. 4 Figure of weighted smoothness of reconstructed data at different scales

图5 为1 至6 尺度组合的阈值重构图，属精细信息图件类别，视觉效应上基本等同于原始数据的精细程度，与原始数据相比，仅仅去除了条带，没有产生其他降噪方法经常出现的边界效应，测区边界数据没有畸变，信息损失较少，可以做局部异常特征分析及异常解释及1∶5 万标准图件。

图5 某测区钍元素含量曲波变换第1 至第6 尺度合并阈值重构图Fig. 5 The map reconstructed by the curve wave transform of scale 1 to 6 combined with threshold values of the content of thorium element in study area

图6 为1 至4 尺度组合的阈值重构图，属细节信息图件类别，视觉效应上去除了条带，碎小斑块较少。可以用来做航空伽马能谱异常分布规律及分类解释，一方面信息足够，另一方面可以不受特别碎小斑块的影响。

图6 某测区钍元素含量曲波处理第1 至第4 尺度合并阈值重构图Fig. 6 The map reconstructed by the curve wave transform of scale 1 to 4 combined with threshold values of the content of thorium element in study area

图7 为1 至2 尺度组合的阈值重构图，属粗糙信息图件类别，视觉效应上去除了条带，特征场信息特别清晰，可用来做区域场特征分析、迁移富集规律分析、补径排特征分析、报告插图等。

图7 某测区钍元素含量曲波处理第1 和第2 尺度合并阈值重构图Fig. 7 The map reconstructed by the curve wave transform of scale 1 and 2 combined with threshold values of the content of thorium element in study area

综上所述，所有处理结果均消除了沿测线方向不同尺度（315 度方向）的条带，降低了随机噪声，没有产生边界效应，测区边界数据没有畸变，具有较好的效果。随着所用尺度的增加，反应的信息也越来越精细，不同尺度的阈值重构图在数据解释方面可以发挥不同的作用。

3 结 论

综上所述，曲波变换用于航空伽马能谱数据处理具有如下优势：

1）利用曲波变换多尺度的特点，可以将航空伽马能谱重构为不同尺度的图件，适用于不同的数据解释。粗尺度图件适用于区域场特征分析、迁移富集规律分析、补径排特征分析等；细节尺度图件适用于航空伽马能谱异常分布规律及分类解释；精细尺度图件适用于局部异常特征分析及异常解释。

2）利用曲波变换多尺度和多方向的特点，可以快捷方便地去除其他降噪方法非常难于去除的沿测线方向不同尺度的条带，并且不降低分辨率，不损失局部异常信息。

3）在降低随机噪声时，不产生边界效应，测区边界数据没有畸变。

实验结果表明：曲波变换用于航空伽马能谱数据处理快捷方便、信息损失少、成果丰富，效果较好。