CN117607804B - 一种基于f-k变换改进的探地雷达杂波抑制方法 - Google Patents
一种基于f-k变换改进的探地雷达杂波抑制方法 Download PDFInfo
- Publication number
- CN117607804B CN117607804B CN202410093334.1A CN202410093334A CN117607804B CN 117607804 B CN117607804 B CN 117607804B CN 202410093334 A CN202410093334 A CN 202410093334A CN 117607804 B CN117607804 B CN 117607804B
- Authority
- CN
- China
- Prior art keywords
- matrix
- interference
- data
- domain
- clutter
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 230000000149 penetrating effect Effects 0.000 title claims abstract description 42
- 230000009466 transformation Effects 0.000 title claims abstract description 17
- 230000001629 suppression Effects 0.000 title claims abstract description 14
- 230000006872 improvement Effects 0.000 title claims description 8
- 230000008878 coupling Effects 0.000 claims abstract description 8
- 238000010168 coupling process Methods 0.000 claims abstract description 8
- 238000005859 coupling reaction Methods 0.000 claims abstract description 8
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 8
- 238000004458 analytical method Methods 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 61
- 230000000694 effects Effects 0.000 claims description 17
- 230000008569 process Effects 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 7
- 239000012535 impurity Substances 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 4
- 230000001131 transforming effect Effects 0.000 claims description 4
- 238000001228 spectrum Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 2
- 238000007781 pre-processing Methods 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims description 2
- 238000001914 filtration Methods 0.000 abstract description 3
- 230000003044 adaptive effect Effects 0.000 abstract description 2
- 238000006243 chemical reaction Methods 0.000 abstract 1
- 239000002689 soil Substances 0.000 description 6
- 238000010586 diagram Methods 0.000 description 4
- 238000013135 deep learning Methods 0.000 description 3
- 230000000877 morphologic effect Effects 0.000 description 3
- 238000013528 artificial neural network Methods 0.000 description 2
- 230000002146 bilateral effect Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000002592 echocardiography Methods 0.000 description 2
- 238000010801 machine learning Methods 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000012549 training Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000009933 burial Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002372 labelling Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/023—Interference mitigation, e.g. reducing or avoiding non-intentional interference with other HF-transmitters, base station transmitters for mobile communication or other radar systems, e.g. using electro-magnetic interference [EMI] reduction techniques
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/36—Means for anti-jamming, e.g. ECCM, i.e. electronic counter-counter measures
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明提出了一种基于F‑K变换改进的探地雷达杂波抑制方法。首先，通过二维傅里叶变换，将探地雷达B‑scan时域数据转换到F‑K域上进行分析。接着对F‑K域数据进行中心化，直达波以及地表反射强信号集中在频率轴附近。使用自适应的SVD奇异值分解方法，在F‑K域中滤除直耦波和地表干扰。经过实际杂波情况，计算并设计出极坐标下边缘高斯缓变的F‑K滤波器，滤除线性干扰。最后将滤波后的数据，去中心化，并且二维傅里叶反变换，回到时域上，即可得到抑制杂波后的B‑scan结果。
Description
技术领域
本发明涉及雷达领域，尤其涉及一种基于F-K变换改进的探地雷达杂波抑制方法。
背景技术
探地雷达（Ground penetrating radar,GPR）是超宽带雷达的一种，其工作原理是通过向地下或目标媒质发射电磁波，利用电磁波在不同电性的介质中传播规律不同的特点，通过对接收到的雷达回波进行分析得到目标媒质的电磁参数和结构，从而实现对介质内部掩埋目标的非侵入式探测。
得益于20世纪80年代和20世纪90年代计算机技术的发展，探地雷达信号的采集与处理效率得到提高。Birken研发出支持多通道数据采集的探地雷达系统；Chew开发出用于仿真3维探地雷达回波数据的软件；Nuzzo实现3维探地雷达数据的可视化；Roberts分析了全极化探地雷达系统的极化特征。
在实际测试的过程中，由于探地雷达系统辐射出的电磁波需要穿透空气-土壤分界面以及地下土壤介质才能到达目标位置处，因此来自地下感兴趣目标的回波通常淹没在杂波下。如何抑制探地雷达回波数据中的杂波信号，一直是探地雷达信号处理领域的一个研究热点。
近年来，国内外学者提出很多探地雷达数据杂波抑制的信号处理方法，可以简单地分为滤波器方法、杂波建模方法、子空间方法、变换域方法、稀疏表示的方法、机器学习方法。
滤波器方法包括时间选通法（Time-gating）、平均减法（Mean Subtraction，MS）、水平高通滤波器、中值滤波器、数字高通滤波器、对称滤波器、多尺度双边滤波器(Multiscale Directional Bilateral Filter，MDBF)以及格型滤波器。其受限于需要设计特殊的滤波器，设置特定的参数才能进行滤波。且面对复杂的场景难以获得较好的滤波效果。
子空间方法将探地雷达B-scan数据分解为杂波信号子空间和目标信号子空间。由于杂波信号要比目标信号强得多，因此杂波信号子空间在B-scan数据中占主要成分。在子空间中去除杂波成分，留下的便是目标的波形了。
在探地雷达时域回波数据中，杂波信号和目标信号可能无法分离或者难以分离，因此变换域方法的核心思路是将探地雷达时域回波数据变换到其他可以将杂波信号和目标信号相分离的域中。常见的有一维小波变换、二维小波变换、F-K（频率-波数）变换等。
基于稀疏表示的形态成分分析（Morphological Component Analysis，MCA）方法可以产生两个不同形态分量的可相互区分的稀疏表示，这对杂波信号和目标信号的进一步有效的分离是非常有利的。但是其构建稀疏字典的过程较为耗时，最后的杂波抑制效果泛化性较差。
近年兴起的深度学习方法已经在计算机视觉领域取得惊人的效果。与传统机器学习方法不同的是深度学习方法采用神经网络提取探地雷达B-scan数据的高维抽象特征，不需要繁琐的处理流程。神经网络通过在人工标注B-scan数据中的目标信号来学习区分杂波信号和目标信号，因此需要大量有标注的训练数据，实际应用过程中的效果较为一般。深度学习方法的性能受训练数据集的影响很大，可能存在较差的泛化问题。
但是以上的方法，大多只能针对地表直耦波，或是地面的杂波进行处理，对于复杂环境，综合场景下产生的杂波，难以做到针对性地滤除。尤其是对于线性的干扰，没有较好的抑制方法。
发明内容
基于以上的观点，提出了一种基于F-K变换改进的探地雷达杂波抑制方法。首先，通过二维傅里叶变换，将探地雷达B-scan时域数据转换到F-K域上进行分析。接着对F-K域数据进行中心化，直达波以及地表反射强信号集中在频率轴附近。使用自适应的SVD奇异值分解方法，在F-K域中滤除直耦波和地表干扰。经过实际杂波情况，计算并设计出极坐标下边缘高斯缓变的F-K滤波器，滤除线性干扰。最后将滤波后的数据，去中心化，并且二维傅里叶反变换，回到时域上，即可得到抑制杂波后的B-scan结果。
本发明提出一种基于F-K变换改进的探地雷达杂波抑制方法，所述方法包括以下步骤：
步骤S1，对原始采集到的探地雷达B_scan数据进行去背景预处理，获得矩阵M；
步骤S2，对矩阵M进行二维傅里叶变换，变换到F-K域上，进行中心化频率偏移的处理，获得矩阵N；
步骤S3，对矩阵N进行SVD奇异值分解，得到矩阵N的奇异值；
步骤S4，对奇异值使用一维k-means的方法，进行聚类的处理聚成n类，去除直达波分量干扰以及抖动分量干扰后，经过奇异值反变换相乘后得到去除强干扰后的数据/>；
步骤S5，根据实际应用场景，设计位于F-K域下，极坐标系中的边缘高斯缓变的滤波器；
步骤S6，将处理过的矩阵与设计的边缘高斯缓变的滤波器线性相乘，得到去除了线性干扰后的数据/>；
步骤S7，将去除杂波干扰后的矩阵去中心化，通过二维的傅里叶逆变换，得到去除杂波后的时域B-scan数据/>。
进一步，步骤S1的具体实现方法为，用探地雷达系统采集数据；设原始的B-scan数据由x条A-scan数据组成，每条A-scan数据中有t个样值点，记作矩阵，对收集到的B_scan数据在时域上进行预处理，使用慢时间去背景的方法，去除直耦的杂波、固定的背景以及多次反射的干扰，
处理过程如下：
(1)
上式通过在慢时间方向上减去均值，达到去除固定背景的效果，其中，为原始采样的B_scan数据，M为预处理后的矩阵，i为第i条A-scan数据，其取值范围为1到x，j为第j个样值点，其取值范围为1到t。
进一步，步骤S2的具体实现方法为：
将预处理过的矩阵M，通过二维的傅里叶变换，变换到F-K域上；
(2)
其中，变换后F-K域上的B-scan数据即为，f表示频率，k表示波数，为了后续处理更加的方便，在矩阵/>的中心位置建立坐标系，将/>的第一象限与第三象限交换，将第二象限与第四象限交换，以达成将零频分量移到频谱中心的效果，频率偏移过后的矩阵记作N。
进一步，步骤S3的具体实现方法为，对中心化后的矩阵N进行SVD奇异值分解：
(3)
其中，是f阶的酉矩阵，/>是k阶的酉矩阵，/>，矩阵S仅主对角线上有值，为矩阵N的奇异值/>；
上式也可以写成：
(4)
其中，为矩阵/>的第i个列向量，/>为矩阵/>的第j个行向量，/>为矩阵N的第k个奇异值；
由探地雷达原理进行分析，矩阵N分为至少四个分量：
(5)
其中，代表了地面或是墙面的一次反射波分量，/>代表了地面起伏造成的杂波分量，/>代表了线性的杂波干扰分量，/>代表了目标的分量。
进一步，步骤S4的具体实现方法为：
取出S矩阵中对角的奇异值，将其排列为一维的数组：
(6)
对奇异值数组使用一维k-means的方法，进行聚类的处理聚成n类，具体的方法流程如下所述：
步骤S41：进行中心点初始化，为n个聚类分组选择初始的中心点；
步骤S42：计算这r个对象与n个中心点的欧式距离；
步骤S43：对比每个对象与中心点的距离，将每个对象分配给距离最近的中心点，从而达到聚类的效果；
步骤S44：根据聚类分组中的样本数值，重新计算平均值，更新每个聚类的中心点；
步骤S45：重复步骤S42、S43、S44，直到聚类算法收敛，即中心点和聚类的分组多次迭代后不再发生改变，停止上述的迭代；
最后能将聚为四类其分别对应N的各个分量：
(7)
由于奇异值分量较大的两类，分别为a个地面墙面的直达分量以及b个地面抖动的强干扰分量，将前两个分量的奇异值直接清零，得到了：
(8)
将奇异值对角矩阵S中的主对角线值用进行替换，最后则得到了去除直达波分量干扰以及抖动干扰后的/>，经过奇异值反变换相乘后得到去除强干扰后的数据/>：
(9)
进一步，步骤S5的具体实现方法为：
设计合适的F-K滤波器；首先，视速度在时域的定义如下：
(10)
其中，是入射电磁波的速度，/>是电磁波的入射角；
在F-K域中，数据点处的视速度/>定义为/>，即为/>数据点与F-K域圆心/>连线的斜率；
经过模型的计算，得到干扰目标的视速度范围，将直角坐标的斜率映射到极坐标的角度，得到杂波的角度范围；
采用边缘高斯缓变的滤波器：
(11)
其中，的取值范围为/>，根据测试经验值，将标准差/>设置为0.1，系数/>取/>，均值/>设置为/>。
进一步，步骤S7的具体实现方法为，将去除杂波干扰后的矩阵第一象限与第三象限交换，将第二象限与第四象限交换，达到去中心化的效果，f表示频率，k表示波数，通过二维的傅里叶逆变换，回到时域的B-scan数据/>，x为A-scan条数，t为样值点数：
(12)
即为最后处理得到的，去除杂波后的数据结果。
本发明具有以下有益技术效果：
（1）本发明提出了一种探地雷达信号杂波抑制方法，能够有效去除地表直达波干扰，地表起伏造成的杂波，同时可以有效抑制与目标视速度不同的线性干扰波形。引入了高斯缓变的滤波器边缘，在有效去杂的同时减少了振铃效应导致的图像失真。经过本发明的处理，去除了绝大部分的杂波，使研究人员能更清晰地对探地雷达B-scan数据进行判读。
（2）本发明使用的场景较为的广泛，不仅能够应用在地耦探地雷达上，还可以应用在空耦探地雷达上。针对不同的应用场景，均能有效对探地雷达杂波信号进行有效去除。
（3）本发明模型均为物理模型计算得到的，最后的可解释性较强。可以对探地雷达B-scan有效地处理，使得目标杂波比提升，更容易获取到目标波形。
（4）本发明提出的算法，算法复杂度较低，运算处理的速度较快。所需的算力极少，可部署性较强，应用场景较广。
附图说明
图 1 为本发明的方法流程图；
图 2 为地耦实际示意图；
图 3 为空耦实际示意图；
图 4 为极坐标下边缘高斯缓变的F-K滤波器示意图。
具体实施方式
为了进一步阐明本发明的技术方案、实施要点以及所具备的优势，以下详细说明本发明的实施步骤。在本发明中，提出了一种基于F-K变换改进的探地雷达杂波抑制方法。图1给出了算法流程图。
上述方法的具体实现过程包括：
步骤S1：用所设计的探地雷达系统采集数据，设原始的B-scan数据由x条A-scan数据组成，每条A-scan数据中有t个样值点，记作矩阵。对收集到的B-scan数据在时域上进行预处理。使用慢时间去背景的方法，去除直耦的杂波、固定的背景以及多次反射的干扰。
处理过程如下：
(1)
上式通过在慢时间方向上减去均值，达到去除固定背景的效果。预处理过后的矩阵记为M。其中，为原始采样的B_scan数据，M为预处理后的矩阵，i为第i条A-scan数据，其取值范围为1到x，j为第j个样值点，其取值范围为1到t。
步骤S2：将预处理过的B-scan时域数据M，通过二维的傅里叶变换，变换到F-K（频率-波数）域上。
(2)
其中，变换后F-K域上的B-scan数据即为。f表示频率，k表示波数，为了后续处理更加的方便，在矩阵/>的中心位置建立坐标系，将/>的第一象限与第三象限交换，将第二象限与第四象限交换，以达成将零频分量移到频谱中心的效果。频率偏移过后的矩阵记作N。
步骤S3：对中心化后的矩阵N进行SVD奇异值分解：
(3)
其中，是f阶的酉矩阵，/>，/>是k阶的酉矩阵，矩阵S仅主对角线上有值，为矩阵N的奇异值/>。
上式也可以写成：
(4)
其中，为矩阵/>的第i个列向量，/>为矩阵/>的第j个行向量，/>为矩阵N的第k个奇异值；
由探地雷达原理进行分析，矩阵N可以分为至少四个分量：
(5)
其中，代表了地面或是墙面的一次反射波分量，/>代表了地面起伏造成的杂波分量，/>代表了线性的杂波干扰分量，/>代表了目标的分量。由目标的分量经过介质的衰减，最后的幅度较小，为弱信号。地面或墙面一次反射分量以及地面起伏造成的杂波分量未经过土壤或墙面内介质的衰减，为强信号。所以/>和/>占据了奇异值分解中，奇异值较大的分量。
步骤S4：取出S矩阵中对角的奇异值，将其排列为一维的数组：
(6)
对奇异值数组使用一维k-means的方法，进行聚类的处理聚成n类，具体的方法流程如下所述：
步骤S41：进行中心点初始化，为n个聚类分组选择初始的中心点，这些中心点称为Means，可以依据经验，也可以用随机数随意选择；
步骤S42：计算这r个对象与n个中心点的欧式距离；
步骤S43：对比每个对象与中心点的距离，将每个对象分配给距离最近的中心点，从而达到聚类的效果；
步骤S44：根据聚类分组中的样本数值，重新计算平均值，更新每个聚类的中点。
步骤S45：重复步骤S42、S43、S44，直到聚类算法收敛，即中心点和聚类的分组多次迭代后不再发生改变，停止上述的迭代。
最后能将聚为四类其分别对应N的各个分量：
(7)
由于奇异值分量较大的两类，分别为a个地面或墙面的直达分量以及b个地面抖动的强干扰分量，所以将前两个分量的奇异值直接清零。得到了：
(8)
将奇异值对角矩阵S中的主对角线值用进行替换，最后则得到了去除直达波分量干扰以及抖动干扰后的/>。经过奇异值反变换相乘后得到去除强干扰后的数据/>：
(9)
步骤S5：经过理论推导，设计合适的F-K滤波器。首先，视速度在时域的定义如下：
(10)
其中，是入射电磁波的速度，/>是电磁波的入射角。
在F-K域中，数据点处的视速/>定义为/>，即为/>数据点与F-K域圆心/>连线的斜率。
针对地耦探地雷达的实际情景，如图2所示：
图中T为探地雷达的发射天线，R为接收天线。对于埋深不同的干扰目标A和待测目标B，最后可以得到不同的探测角度，/>，继而推出不同的视速度/>和/>。这个是在F-K域上分离干扰和目标的重要依据。
针对空耦探地雷达的实际情况，如图3所示：
由上述地耦的推导推广到空耦雷达的适用场景，空耦雷达距离地面h的高度，假设空气中的介电常数为，土层的相对介电常数为/>。一空气中的目标，距离雷达h1，其夹角为/>。土地内有一目标深度为h2，其体现出来的视速角度为/>。由于目标在土层之中，会产生折射，其中，/>为在介质中折射后，实际的视速角度。
经过模型的计算，可以得到干扰目标的视速度范围，将直角坐标的斜率映射到极坐标的角度，可以得到杂波的角度范围。
为了防止直接设计滤波器，边缘跳变严重，导致最后的B-scan数据有严重的振铃现象，本发明设计了边缘高斯缓变的滤波器。
(11)
其中，的取值范围为/>，根据测试经验值，将标准差/>设置为0.1，系数/>取/>，均值/>设置为/>。
图4为一个极坐标下边缘高斯缓变的F-K滤波器示意图。
步骤S6：将处理过的矩阵与设计的滤波器线性相乘，得到去除了线性干扰后的数据/>。
步骤S7：将去除杂波干扰后的矩阵第一象限与第三象限交换，将第二象限与第四象限交换，达到去中心化的效果。通过二维的傅里叶逆变换，回到时域的B-scan数据：
(12)
其中，x为A-scan条数，t为样值点数：即为最后处理得到的，去除杂波后的数据结果。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (7)
1.一种基于F-K变换改进的探地雷达杂波抑制方法，其特征在于，所述方法包括以下步骤，
步骤S1，对原始采集到的探地雷达B_scan数据进行去背景预处理，获得矩阵M；
步骤S2，对矩阵M进行二维傅里叶变换，变换到F-K域上，进行中心化频率偏移的处理，获得矩阵N；
步骤S3，对矩阵N进行SVD奇异值分解，得到矩阵N的奇异值；
步骤S4，对奇异值使用一维k-means的方法，进行聚类的处理聚成n类，去除直达波分量干扰以及抖动分量干扰后，经过奇异值反变换相乘后得到去除强干扰后的数据/>；
步骤S5，根据实际应用场景，设计位于F-K域下，极坐标系中的边缘高斯缓变的滤波器；
步骤S6，将处理过的矩阵与设计的边缘高斯缓变的滤波器线性相乘，得到去除了线性干扰后的数据/>；
步骤S7，将去除杂波干扰后的矩阵去中心化，通过二维的傅里叶逆变换，得到去除杂波后的时域B-scan数据/>。
2.根据权利要求1所述的方法，其特征在于，步骤S1的具体实现方法为，用探地雷达系统采集数据；设原始的B-scan数据由x条A-scan数据组成，每条A-scan数据中有t个样值点，记作矩阵，对收集到的B_scan数据在时域上进行预处理，使用慢时间去背景的方法，去除直耦的杂波、固定的背景以及多次反射的干扰；
处理过程如下：
(1)
上式通过在慢时间方向上减去均值，达到去除固定背景的效果，其中，为原始采样的B_scan数据，M为预处理后的矩阵，i为第i条A-scan数据，其取值范围为1到x，j为第j个样值点，其取值范围为1到t。
3.根据权利要求2所述的方法，其特征在于，步骤S2的具体实现方法为：
将预处理过的矩阵M，通过二维的傅里叶变换，变换到F-K域上：
(2)
其中，变换后F-K域上的B-scan数据即为，f表示频率，k表示波数，为了后续处理更加的方便，在矩阵/>的中心位置建立坐标系，将/>的第一象限与第三象限交换，将第二象限与第四象限交换，以达成将零频分量移到频谱中心的效果，频率偏移过后的矩阵记作N。
4.根据权利要求3所述的方法，其特征在于，步骤S3的具体实现方法为，对中心化后的矩阵N进行SVD奇异值分解：
(3)
其中，是f阶的酉矩阵，/>是k阶的酉矩阵，/>，矩阵S仅主对角线上有值，为矩阵N的奇异值/>；
上式也可以写成：
(4)
其中，为矩阵/>的第i个列向量，/>为矩阵/>的第j个行向量，/>为矩阵N的第k个奇异值；
由探地雷达原理进行分析，矩阵N分为至少四个分量：
(5)
其中，代表了地面或是墙面的一次反射波分量，/>代表了地面起伏造成的杂波分量，/>代表了线性的杂波干扰分量，/>代表了目标的分量。
5.根据权利要求4所述的方法，其特征在于，步骤S4的具体实现方法为：
取出S矩阵中对角的奇异值，将其排列为一维的数组：
(6)
对奇异值数组使用一维k-means的方法，进行聚类的处理聚成n类，具体的方法流程如下所述：
步骤S41：进行中心点初始化，为n个聚类分组选择初始的中心点；
步骤S42：计算这r个对象与n个中心点的欧式距离；
步骤S43：对比每个对象与中心点的距离，将每个对象分配给距离最近的中心点，从而达到聚类的效果；
步骤S44：根据聚类分组中的样本数值，重新计算平均值，更新每个聚类的中心点；
步骤S45：重复步骤S42、S43、S44，直到聚类算法收敛，即中心点和聚类的分组多次迭代后不再发生改变，停止上述的迭代；
最后能将聚为四类其分别对应N的各个分量：
(7)
由于奇异值分量较大的两类，分别为a个地面墙面的直达分量以及b个地面抖动的强干扰分量，将前两个分量的奇异值直接清零，得到了：
(8)
将奇异值对角矩阵S中的主对角线值用进行替换，最后则得到了去除直达波分量干扰以及抖动干扰后的/>，经过奇异值反变换相乘后得到去除强干扰后的数据/>：
(9)。
6.根据权利要求5所述的方法，其特征在于，步骤S5的具体实现方法为：
设计合适的F-K滤波器；首先，视速度在时域的定义如下：
(10)
其中，是入射电磁波的速度，/>是电磁波的入射角；
在F-K域中，数据点处的视速度/>定义为/>，即为/>数据点与F-K域圆心/>连线的斜率；
经过模型的计算，得到干扰目标的视速度范围，将直角坐标的斜率映射到极坐标的角度，得到杂波的角度范围；
采用边缘高斯缓变的滤波器：
(11)
其中，的取值范围为/>，根据测试经验值，将标准差/>设置为0.1，系数/>取，均值/>设置为/>。
7.根据权利要求6所述的方法，其特征在于，步骤S7的具体实现方法为，将去除杂波干扰后的矩阵第一象限与第三象限交换，将第二象限与第四象限交换，达到去中心化的效果，f表示频率，k表示波数，通过二维的傅里叶逆变换，回到时域的B-scan数据/>，x为A-scan条数，t为样值点数：
(12)
即为最后处理得到的，去除杂波后的数据结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410093334.1A CN117607804B (zh) | 2024-01-23 | 2024-01-23 | 一种基于f-k变换改进的探地雷达杂波抑制方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410093334.1A CN117607804B (zh) | 2024-01-23 | 2024-01-23 | 一种基于f-k变换改进的探地雷达杂波抑制方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117607804A CN117607804A (zh) | 2024-02-27 |
| CN117607804B true CN117607804B (zh) | 2024-03-22 |
Family
ID=89956563
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410093334.1A Active CN117607804B (zh) | 2024-01-23 | 2024-01-23 | 一种基于f-k变换改进的探地雷达杂波抑制方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117607804B (zh) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2440590A1 (en) * | 2002-09-12 | 2004-03-12 | Kelman Technologies Inc. | Method of using matrix rank reduction to remove random noise from seismic data |
| CN109117784A (zh) * | 2018-08-08 | 2019-01-01 | 上海海事大学 | 一种改进经验模态分解的船舶电力推进系统故障诊断方法 |
| CN113189641A (zh) * | 2021-03-25 | 2021-07-30 | 西安石油大学 | 一种两道多模式瑞利波地下探测系统及方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6895336B2 (en) * | 2002-09-12 | 2005-05-17 | Kelman Technologies, Inc. | Method of using matrix rank reduction to remove random noise from seismic data |
| CN112505648B (zh) * | 2020-11-19 | 2023-06-30 | 西安电子科技大学 | 基于毫米波雷达回波的目标特征提取方法 |
-
2024
- 2024-01-23 CN CN202410093334.1A patent/CN117607804B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2440590A1 (en) * | 2002-09-12 | 2004-03-12 | Kelman Technologies Inc. | Method of using matrix rank reduction to remove random noise from seismic data |
| CN109117784A (zh) * | 2018-08-08 | 2019-01-01 | 上海海事大学 | 一种改进经验模态分解的船舶电力推进系统故障诊断方法 |
| CN113189641A (zh) * | 2021-03-25 | 2021-07-30 | 西安石油大学 | 一种两道多模式瑞利波地下探测系统及方法 |
Non-Patent Citations (3)
| Title |
|---|
| f-k滤波在探地雷达数据去噪中的应用;吴海波;张平松;胡雄武;秦镇;;安徽理工大学学报(自然科学版);20200115(第01期);全文 * |
| 基于探地雷达的地下管线数据处理研究;张力文;中国优秀硕士学位论文全文数据库 基础科学辑 (月刊);20220515(第05期);全文 * |
| 浅谈地质雷达数据的精细处理;苏智光;钱东宏;廖建军;李仕胜;;工程勘察;20120701(第07期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117607804A (zh) | 2024-02-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108872962B (zh) | 基于分数阶傅里叶变换的激光雷达微弱信号提取和分解方法 | |
| CN109581516B (zh) | 曲波域统计量自适应阈值探地雷达数据去噪方法及系统 | |
| CN110133643B (zh) | 植物根系探测方法及装置 | |
| CN108985304B (zh) | 一种基于浅剖数据的沉积层结构自动提取方法 | |
| CN109669182B (zh) | 无源双基地sar动/静目标联合稀疏成像方法 | |
| CN109143363B (zh) | 海洋拖缆双检采集鬼波压制方法及系统 | |
| CN108710888B (zh) | 一种探地雷达数据配准方法 | |
| CN108983287B (zh) | 一种基于凸集投影算法的曲波变换反假频地震数据重建方法 | |
| CN104077749A (zh) | 一种基于轮廓波变换的地震数据去噪方法 | |
| CN108802725A (zh) | 一种浅层穿透雷达合成孔径成像方法 | |
| CN112255607A (zh) | 一种海杂波的抑制方法 | |
| CN112327362A (zh) | 速度域的海底多次波预测与追踪衰减方法 | |
| CN109581481B (zh) | 一种便携式高频可控震源地震信号谐波干扰消除方法 | |
| CN117607804B (zh) | 一种基于f-k变换改进的探地雷达杂波抑制方法 | |
| Zhao et al. | Seismic data denoising using curvelet transforms and fast non-local means | |
| CN108919345B (zh) | 一种海底电缆陆检噪声的衰减方法 | |
| CN109975873B (zh) | 一种逆时偏移成像去除低频噪音的方法及系统 | |
| CN116203634A (zh) | 一种基于低秩约束的鬼波去除方法 | |
| CN110703332B (zh) | 一种鬼波压制方法 | |
| CN116224324A (zh) | 基于深度学习的超分辨率3d-gpr图像的频率-波数分析方法 | |
| CN109427042B (zh) | 一种提取局部海域沉积层的分层结构和空间分布的方法 | |
| CN111257938A (zh) | 基于小波互相关时移地震虚拟震源波场重构方法和系统 | |
| CN116047504A (zh) | 一种改进反褶积抑制探地雷达多次波的方法 | |
| CN115113163A (zh) | 一种探地雷达多分辨率低秩稀疏分解杂波抑制方法 | |
| Li et al. | Multiple attention mechanisms-based convolutional neural network for desert seismic denoising |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |