CN107450054A - 一种自适应探地雷达数据去噪方法 - Google Patents

Abstract

本发明公开了一种自适应探地雷达数据去噪方法，通过将探地雷达数据利用奇异值分解成一系列与不同奇异值对应的数据成分，以便于根据不同信号或噪声类型单独处理数据；然后利用探地雷达数据的初至同相轴近似探测对象表面起伏状况，同时计算初至同相轴的均方根高；最后将均方根高作为自动判断与相干噪声对应的奇异值的依据，从而达到自动去除相干噪声的目的，避免了基于SVD分解方法需要人为判断奇异值的麻烦，直接根据均方根高参数来判断奇异值，既准确又省时。

Description

一种自适应探地雷达数据去噪方法

技术领域

本发明属于工程勘探数据处理技术领域，具体涉及一种自适应探地雷达数据去噪方法。

背景技术

探地雷达自上世纪初Letmbach和Lowy提出概念至今已走过1个多世纪的历程，不仅理论成熟，实际应用也取得了巨大成果。然而在实际应用中，雷达信号不可避免地会受到各种噪声干扰，从而导致信噪比降低，深部信号被掩盖，成像效果差，由雷达信号反演计算各类参数的精确性更是会受到不同程度的影响，使得最终成果质量大幅降低。

从噪声来源来分，探地雷达噪声主要可以分为仪器噪声、天线耦合噪声、人文工业噪声、介质内的非均匀散射、强地表反射、起伏地表散射以及各种多次反射波等；针对各类噪声，已有的去噪方法主要包括似然比检测技术、参数***识别技术、小波包分解技术、子空间滤波技术、独立成分分析、平均值滤波技术、频率域滤波技术等，虽然这些方法有充分的理论依据，在模型数据上也能获得理想效果，但受限于各种假设和实际条件，在处理受噪声污染的实际数据时通常难以达到预期目标。

因此，寻求高精度、高效率、适应性强的去噪方法是成功应用探地雷达解决实际问题的前提基础。奇异值分解(SVD)是以线性代数理论为基础的矩阵分解方法，其核心思想是将高维复杂矩阵表示成低维简单矩阵，从而突出原矩阵的主要特征；SVD技术在各个领域均有广泛应用，包括图像压缩、信号去噪、模式识别、特征提取等，然而这些应用都面临一个共同的问题，即如何准确决定奇异值个数，而这是成功利用SVD技术的关键步骤。为此，各个领域的专家学者根据特定目标发展了各自的处理方法，也取得了一定成效，然而还没有哪种方法具备广泛的适应性。

在探地雷达去除水平规则干扰方面，Abujarad等简单地将第一个奇异值所对应成分视为水平干扰，将剩余数据均作为有效信号，该方法只对极其规则的水平干扰噪声有效。此外，还包括奇异熵增量、信噪比经验等方法，这些方法都只是半定量方法，通常依赖于处理人员的经验，这不仅降低了工作效率，而且缺乏客观性。针对此问题，亟需在SVD分解的基础上研究自动判断奇异值个数的方法，并根据奇异值重新合成目标成分来实现准确高效的去噪目的，从而提高雷达数据的信噪比，突出目标信号特征，为后续处理提供高质量的基础数据，保证结果的可靠性。

发明内容

鉴于上述，本发明提供了一种自适应探地雷达数据去噪方法，能够自动确定相干噪声奇异值个数，避免针对不同数据需人机交互判定奇异值数目带来的麻烦。

一种自适应探地雷达数据去噪方法，包括如下步骤：

(1)对于原始探地雷达数据，截取其前k纳秒内的数据并确定数据的标准道，k为大于0的实数；

(2)对于任一扫描道，计算该扫描道与标准道在前k纳秒内的互相关函数并取其最大值；

(3)根据步骤(2)遍历所有扫描道，得到由各扫描道互相关函数最大值组成的初至曲线L，进而计算出初至曲线L的均方根高D；

(4)根据以下公式计算出与相干噪声对应的奇异值个数N；

N＝ROUND{0.2634D+1.3086}

其中：ROUND{}为向上取整函数；

(5)对原始探地雷达数据进行奇异值分解，剔除奇异值分解后的前N个数据成分，将剩余数据成分重新合成后即完成去噪过程。

所述步骤(1)中的原始探地雷达数据为n×m大小的数据矩阵，n和m分别对应为时间和空间上的采样点数；所述扫描道即对应探测对象表面上任一采样点的雷达数据，所述标准道即对应探测对象表面上某一采样点的雷达数据，该采样点为随机或人为指定。

所述步骤(2)中的互相关函数表达式如下：

其中：R_j(τ)为第j扫描道与标准道在前k纳秒内的互相关函数，x₀(t)为t时刻的标准道，t和τ均为时间采样序号且1≤τ≤K，K为前k纳秒内雷达数据的采样点数，x_j ^*(t-τ)为x_j(t-τ)的共轭复数，x_j(t-τ)为t-τ时刻的第j扫描道，j为自然数且1≤j≤m。

所述步骤(3)中根据以下公式计算初至曲线L的均方根高D：

其中：d(j)为第j扫描道与标准道在前k纳秒内的互相关函数最大值，为初至曲线L的平均值。

优选地，k＝5。

本发明通过将探地雷达数据利用奇异值分解成一系列与不同奇异值对应的数据成分，以便于根据不同信号或噪声类型单独处理数据；然后利用探地雷达数据的初至同相轴近似探测对象表面起伏状况，同时计算初至同相轴的均方根高；最后将均方根高作为自动判断与相干噪声对应的奇异值的依据，从而达到自动去除相干噪声的目的，避免了基于SVD分解方法需要人为判断奇异值的麻烦，直接根据均方根高参数来判断奇异值，既准确又省时。

附图说明

图1为探测物表面平整情况下获取实际含噪声的探地雷达数据示意图。

图2为图1利用互相关方法自动计算得到的初至曲线结果示意图。

图3为图1利用奇异值分解后得到的前50个奇异值分布结果示意图。

图4为图1利用本发明方法提取到的噪声示意图。

图5为图1利用本发明方法去噪后的数据示意图。

图6为探测物表面非平整情况下获取实际含噪声的探地雷达数据示意图。

图7为图6利用互相关方法自动计算得到的初至曲线结果示意图。

图8为图6利用奇异值分解后得到的前50个奇异值分布结果示意图。

图9为图6利用本发明方法提取到的噪声示意图。

图10为图6利用本发明方法去噪后的数据示意图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

本发明利用奇异值分解将原始探地雷达数据分解成不同成分，利用初至曲线的均方根高作为决定相干噪声分布范围的参数，提取噪声以外奇异值分解成分合成新数据，实现自适应去噪目的。

设A是n×m阶原始探地雷达数据，其中n和m分别是时间和空间采样点数，将A进行奇异值分解，如下式所示：

其中：U是n×n阶正交矩阵，其每列为协方差矩阵AA^T的特征向量u_i(1≤i≤n)，称为A的左奇异矩阵；V是m×m阶正交矩阵，其每列为协方差矩阵A^TA的特征向量v_j(1≤j≤m)，称为A的右奇异矩阵；Σ是一个由r(r是A的秩)个元素组成的n×m阶对角阵，其元素为从大到小排列的σ_s(1≤s≤r)，且σ₁≥σ₂≥...≥σ_r≥0；σ_s称为A的奇异值，实际上就是A的特征值的平方根；将上式展开：

式中每一项u_sv_s ^T是一个与A完全同维的矩阵，它代表了A的第s个特征成分；上式表明A可以被分解成一系列不同的特征成分，σ_s表示权因子，它决定了对应特征成分对构成完整A的贡献大小。根据大量实践经验，绝大多数情况下，前10％的奇异值便占据了所有奇异值之和的99％以上，对于某些特殊情况，前1％的奇异值几乎就控制了整个数据或图像的主要特征；换句话说，可以用利用极少量的奇异值来极大地逼近原数据A，即：

上式中U'_n×k、Σ'_k×k、分别是由U、Σ、V^T的前k个成分组成的矩阵，k通常是一个远小于维度n和m的正整数。

探地雷达探测对象表面越平整，首奇异值越突出，去噪时所选取的奇异值个数越小；反之，表面越起伏不平，奇异值越均匀，去噪时所选取的奇异值个数越多，即去噪时奇异值的选择与探测对象表面起伏情况相关。本发明利用互相关法计算初至曲线来近似起伏表面，利用均方根高来定量刻画初至曲线起伏状况，具体过程如下：

(1)选择原始地质雷达数据前5纳秒时间窗内数据，同时选择一道波形明确、未受噪声干扰的数据作为标准道，利用下式计算标准道与其他各道数据的互相关函数：

式中：j是扫描道号，是x_j()的复共轭，x₀()是选中与其他各道数据进行互相关的标准道，R_j(τ)是对应的互相关函数。标准道x₀()可以随机确定，一般情况下大部分道都可以作为标准道。

(2)得到互相关函数后，再找出其中最大值d；根据互相关函数和初至时间标准计算所有道的初至时间，将各道初至时间按顺序绘图形成曲线作为探地雷达数据的初至同相轴。

d＝max(R_j(τ))

d表示了当前道与标准道的偏移量，求出所有道与标准道的互相关函数的最大值，并将他们按顺序排放，就得到了初至同相轴L，其表达式为：

L＝d(j)，1≤j≤m

(3)得到起伏曲线L后，计算均方根高D，它描述了曲线在垂直方向的平均起伏程度，其表达式为：

式中：为L的平均值。

将均方根高D代入下式：

N＝0.2634D+1.3086

(4)计算得到与相干噪声对应的奇异值个数N，再将前面奇异值分解后的前N个数据成分剔除并将剩余数据成分重新合成新数据，从而得到最终的去除噪声的数据。

实施例一

如图1所示为探测物表面平整情况下获取的实际含噪声探地雷达数据，从图1中可以看出大量噪声赋存于剖面中，压制了有效信号的准确表达。此外从噪声的形态课题可以初步判断，与之对应的探测物表面较为平整，不存在严重起伏。图2是在图1的基础上截取的5纳秒时间窗数据，选择图2中黑色垂直细实线所在道为标准道与其他各道进行互相关计算，并在此基础上自动获取所有道初至时间，得到了图2中黑色粗实线所示初至曲线。可以看出，曲线的起伏变化形态与所有道的初至起跳时间吻合得很好。对初至曲线计算均方根高D，并将其代入公式N＝0.2634D+1.3086且取整数部分得N＝1，即对于该数据，水平相干干扰对应的数据成分为第1个奇异值对应成分。将图1所示原始数据进行奇异值分解，并把前30个奇异值绘成如图3所示结果，图3中第1个奇异值即为噪声对应奇异值。从图3中可以看出，第1个奇异值与其他奇异值相差非常大，这也说明了噪声的集中性，只需将第1个值对应成分去除就能达到较好的去噪效果，可见自动判断的奇异值是合理的。图4是将第1个奇异值对应噪声提取出来的结果，对比图1和图4可知，绝大部分噪声都被成功分离。将剩余数据重新合成新数据，也就是去噪后的数据，如图5所示，有效信号得到了很好的回复，水平向噪声被很好地消除了；整个过程自动执行，无需人为决定参数。

实施例二

图6是噪声情况更加严重的一个实际数据，有效信号除了极少区域可分辨以外，大部分区域都被强振幅的噪声覆盖。而且从噪声的形态特征可以判断，探测物表面起伏比较厉害。这给去噪处理带来了很大困难，一般去噪方法都会失效。图7还是在图6基础上截取的前5纳秒时间窗内数据，黑色垂直细实线为标准道，黑色粗实线为根据标准道和互相关自动计算的初至曲线，可见曲线与起伏波动的数据初至也吻合得很好。对初至曲线计算均方根高D，并将其代入公式N＝0.2634D+1.3086且取整数部分得N＝2，即对于该数据，水平相干干扰对应的数据成分为第1、2个奇异值对应成分。图8是将原始数据进行奇异值分解后得到的前30个奇异值分布情况。与实施例一中的奇异值相比，首奇异值有了很大下降，且与后面奇异值的差距大幅减小，这表面原始数据中噪声的分布情况与其他信号的差别在减小，这无疑增加了去噪难度。但较为明显的是前两个奇异值与后续奇异值之间的差异还是较为明显，这也说明了这两个奇异值对应数据成分与后续数据的差异，实际上也就是横向噪声，这也说明了本实施例自动判定奇异值的方法是有效的。图9是提取前两个奇异值对应噪声，与图6所示原始数据相比，大部分噪声也被提取出来了。图10是将剩余奇异值对应成分合成后的数据，有效信号得到了显著加强，原本被覆盖而无法辨识的信号也得到了恢复。此例说明，对于横向不规则的噪声而言，本发明也具有很好的适应性。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。

Claims (5)

1.一种自适应探地雷达数据去噪方法，包括如下步骤：

(1)对于原始探地雷达数据，截取其前k纳秒内的数据并确定数据的标准道，k为大于0的实数；

(2)对于任一扫描道，计算该扫描道与标准道在前k纳秒内的互相关函数并取其最大值；

(3)根据步骤(2)遍历所有扫描道，得到由各扫描道互相关函数最大值组成的初至曲线L，进而计算出初至曲线L的均方根高D；

(4)根据以下公式计算出与相干噪声对应的奇异值个数N；

N＝ROUND{0.2634D+1.3086}

其中：ROUND{}为向上取整函数；

(5)对原始探地雷达数据进行奇异值分解，剔除奇异值分解后的前N个数据成分，将剩余数据成分重新合成后即完成去噪过程。

2.根据权利要求1所述的自适应探地雷达数据去噪方法，其特征在于：所述步骤(1)中的原始探地雷达数据为n×m大小的数据矩阵，n和m分别对应为时间和空间上的采样点数；所述扫描道即对应探测对象表面上任一采样点的雷达数据，所述标准道即对应探测对象表面上某一采样点的雷达数据，该采样点为随机或人为指定。

3.根据权利要求2所述的自适应探地雷达数据去噪方法，其特征在于：所述步骤(2)中的互相关函数表达式如下：

<mrow> <msub> <mi>R</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>t</mi> <mo>=</mo> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </munderover> <msub> <mi>x</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>x</mi> <mi>j</mi> <mo>*</mo> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>-</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> </mrow>

其中：R_j(τ)为第j扫描道与标准道在前k纳秒内的互相关函数，x₀(t)为t时刻的标准道，t和τ均为时间采样序号且1≤τ≤K，K为前k纳秒内雷达数据的采样点数，x_j ^*(t-τ)为x_j(t-τ)的共轭复数，x_j(t-τ)为t-τ时刻的第j扫描道，j为自然数且1≤j≤m。

4.根据权利要求2所述的自适应探地雷达数据去噪方法，其特征在于：所述步骤(3)中根据以下公式计算初至曲线L的均方根高D：

<mrow> <mi>D</mi> <mo>=</mo> <msqrt> <mfrac> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msup> <mrow> <mo>(</mo> <mi>d</mi> <mo>(</mo> <mi>j</mi> <mo>)</mo> <mo>-</mo> <mover> <mi>d</mi> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mi>m</mi> </mfrac> </msqrt> </mrow>

其中：d(j)为第j扫描道与标准道在前k纳秒内的互相关函数最大值，为初至曲线L的平均值。

5.根据权利要求1所述的自适应探地雷达数据去噪方法，其特征在于：k＝5。