CN114154537A - 一种地质雷达信号去噪方法、装置、设备及介质 - Google Patents

一种地质雷达信号去噪方法、装置、设备及介质 Download PDF

Info

Publication number
CN114154537A
CN114154537A CN202111400693.XA CN202111400693A CN114154537A CN 114154537 A CN114154537 A CN 114154537A CN 202111400693 A CN202111400693 A CN 202111400693A CN 114154537 A CN114154537 A CN 114154537A
Authority
CN
China
Prior art keywords
matrix
denoising
geological radar
original
radar signal
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.)
Pending
Application number
CN202111400693.XA
Other languages
English (en)
Inventor
雷洋
江波
田甜
齐法琳
王兆宁
李耀南
刘恒柏
瞿起明
宋国华
贾飞宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Academy of Railway Sciences Corp Ltd CARS
Infrastructure Inspection Institute of CARS
Beijing IMAP Technology Co Ltd
Original Assignee
China Academy of Railway Sciences Corp Ltd CARS
Infrastructure Inspection Institute of CARS
Beijing IMAP Technology Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Academy of Railway Sciences Corp Ltd CARS, Infrastructure Inspection Institute of CARS, Beijing IMAP Technology Co Ltd filed Critical China Academy of Railway Sciences Corp Ltd CARS
Priority to CN202111400693.XA priority Critical patent/CN114154537A/zh
Publication of CN114154537A publication Critical patent/CN114154537A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Signal Processing (AREA)
  • Computing Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Algebra (AREA)
  • Artificial Intelligence (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本文提供了一种地质雷达信号去噪方法、装置、设备及介质，包括：根据地质雷达信号，确定原始矩阵；计算原始矩阵的特征值和特征向量；根据特征值间的斜率比值，从原始矩阵的特征值中筛选出去噪特征值，并根据去噪特征值确定去噪矩阵；根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号，实现了通过根据相邻特征值间的斜率，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵能够将噪声和有效信号对应的去噪矩阵进行划分，并根据所述去噪矩阵和所述原始矩阵实现有效信号的的数据重构，实现了保留原始信号，不需要进行时频域变换，且对于原始信号的时域长度无限制，并且该方法计算简便、参数简单，适用性较强。

Description

一种地质雷达信号去噪方法、装置、设备及介质
技术领域
本发明涉及信号处理技术领域，可用于地质勘探领域，尤其是一种地质雷达信号去噪方法、装置、设备及介质。
背景技术
在地质雷达实际探测过程中，由于仪器本身、人文、自然环境等多方面因素的影响，采集的数据难以避免存在噪声，这些噪声主要为由外界天线、电缆、工程作业等外部复杂环境和地质雷达系统自身引起的随机噪声以及由广播、电视等在地质雷达频段覆盖范围内的通信系统所引起的射频干扰，它们影响了地质雷达数据对地下真实地质体分布情况的反映，通常需要抑制噪声信号，提取所需的目标回波信号，最大程度地还原地下的真实地质情况。
现有去噪方法存在一定的适用条件和局限性，例如采用傅里叶变换将时域信号变换到频域进行去噪的方法，无法处理局部信号，中值滤波法、F-K滤波法和且触发等方法对于平稳信号噪声来说去噪效果良好，而对于频带混叠的含噪数据的去噪效果不理想，因此，需要一种计算简便、参数简单、适用性强又能够较好的抑制地质雷达噪声信号的去噪方法及系统，来提高地质雷达数据的质量。
发明内容
针对现有技术的上述问题，本文的目的在于，提供一种地质雷达信号去噪方法、装置、设备及介质，以解决现有技术中去噪方法存在使用条件局限性，以及运算过程复杂、参数冗余的问题。
为了解决上述技术问题，本文的具体技术方案如下：
一方面，本文提供一种地质雷达信号去噪方法，包括：
根据所述地质雷达信号，确定原始矩阵；
计算所述原始矩阵的特征值、特征向量即特征值间的斜率比值；
根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵；
根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号。
作为本文的一个实施例，所述计算所述原始矩阵的特征值和特征向量，进一步包括：
对所述原始矩阵进行零均值化处理，得到零均值化矩阵；
根据所述零均值化矩阵，确定所述协方差矩阵；
根据所述协方差矩阵，计算所述原始矩阵的特征值和特征向量。
作为本文的一个实施例，所述对所述原始矩阵进行零均值化处理，得到零均值化矩阵，进一步包括：
根据所述原始矩阵中每一列的平均值，构建平均值矩阵；
利用所述原始矩阵减去所述平均值矩阵，得到所述零均值化矩阵。
作为本文的一个实施例，所述根据所述协方差矩阵，计算所述原始矩阵的特征值和特征向量，进一步包括：
根据所述协方差矩阵、单位矩阵和零向量确定所述原始矩阵的特征值；
将所述特征值依次带入至齐次线性方程组中，计算得到每个所述特征值对应的特征向量，其中，所述齐次方程组由所述协方差矩阵、所述单位矩阵和所述零向量组成。
作为本文的一个实施例，计算特征值间的斜率比值，包括：
将所述特征值按照由大到小的顺序进行排序；
将所述特征值的排序序号作为横坐标，将特征值作为纵坐标，建立坐标系；
在所述坐标系下，计算相邻坐标点之间的斜率；
根据所述相邻特征值间的斜率，计算相邻斜率间的斜率比值，将相邻斜率间的斜率比值作为排序后特征值间的斜率比值。
作为本文的一个实施例，所述根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵，进一步包括：
选取最大斜率比值相关的特征值中最大的排序序号；
将特征值排序中位于选取出的排序序号之前的特征值作为所述去噪特征值；
将所述去噪特征值对应的特征向量组合确定所述去噪矩阵。
作为本文的一个实施例，所述根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号，进一步包括：
根据所述原始矩阵中每一列的平均值，构建平均值矩阵；
利用所述原始矩阵减去所述平均值矩阵，得到所述零均值化矩阵；
将所述去噪矩阵和所述零均值化矩阵相乘，得到降维矩阵；
对所述去噪矩阵进行转置处理，得到转置去噪矩阵；
将所述转置去噪矩阵和所述降维矩阵相乘，再与所述平均值矩阵相加，确定所述去噪地质雷达信号。
另一方面，本文还提供一种地质雷达信号去噪装置，包括：
矩阵化单元，用于根据所述地质雷达信号，确定原始矩阵；
分析单元，用于计算所述原始矩阵的特征值、特征向量及特征值间的斜率比值；
筛选单元，用于根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵；
重构单元，用于根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号。
另一方面，本文还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现任意一项所述地质雷达信号去噪方法。
另一方面，本文还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现任意一项所述地质雷达信号去噪方法。
本文通过根据相邻特征值间的斜率，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵能够将噪声和有效信号对应的去噪矩阵进行划分，并根据所述去噪矩阵和所述原始矩阵实现有效信号的的数据重构，实现了保留原始信号，不需要进行时频域变换，且对于原始信号的时域长度无限制，并且该方法计算简便、参数简单，适用性较强。
为让本文的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。
附图说明
为了更清楚地说明本文实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本文的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1示出了本文实施例一种地质雷达信号去噪的整体系统图；
图2示出了本文实施例一种地质雷达信号去噪方法的步骤示意图；
图3示出了本文实施例一种地质雷达信号去噪装置的示意图；
图4示出了本文实施例一种地质雷达信号去噪装置的重构单元示意图；
图5示出了本文实施例一种地质雷达信号去噪方法的数据流程图；
图6示出了本文实施例一种地质雷达信号去噪方法的相邻特征值间的斜率比值示意图；
图7示出了本文实施例一种地质雷达信号去噪方法的筛选条件示意图；
图8示出了本文实施例一种计算机设备示意图。
附图符号说明：
101、地质雷达；
102、运算服务器；
103、控制终端；
301、矩阵化单元；
302、分析单元；
303、筛选单元；
304、重构单元；
3041、平均值模块；
3042、零均值化模块；
3043、降维模块；
3044、转置模块；
3045、运算模块；
802、计算机设备；
804、处理器；
806、存储器；
808、驱动机构；
810、输入/输出模块；
812、输入设备；
814、输出设备；
816、呈现设备；
818、图形用户接口；
820、网络接口；
822、通信链路；
824、通信总线。
具体实施方式
下面将结合本文实施例中的附图，对本文实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本文一部分实施例，而不是全部的实施例。基于本文中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文保护的范围。
需要说明的是，本文的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本文的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、装置、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
如图1所示一种地质雷达信号去噪的整体系统图，包括地质雷达101、运算服务器102和控制终端103。
地质雷达101可以获取地层中的信号，例如可以获取人为爆破产生的地震波，也可以接收施工人员使用声波发送器进行声波发送，并返回来的地质雷达信号，在仪器使用过程中，由于仪器本身或者是使用不规范，以及外界干扰，该外界干扰可以是手机电磁波、无线电发射塔，或者电缆存在电磁干扰等在地质雷达频段覆盖范围内的通信系统所引起的射频干扰，这些因素影响了地质雷达101接收的地质雷达信号的准确性，进行影响了地质雷达信号对应的数据对于地下真实地质体分布情况的反映。
运算服务器102接收地质雷达101的地质雷达信号，可以将声波形式的地质雷达信号进行数字化或者矩阵化，通过这种方式，可以令用户更直观的感受到地下的真实地质情况。在运算服务器102中可以存在傅里叶变化算法，例如FFT(快速傅里叶变换)，通过这种方法，可以实现将地质雷达信号由时域向频域的转化，便于研究地质雷达信号中的功率谱，使得用户通过能量去研究信号，使得研究变得更加简单，令计算高效，但是在使用傅里叶变换时，需要提取利用地质雷达信号的全部时域信息，这是一种整体变换，缺少时域中的定位功能，因此用户需要根据需要使用运算服务器102中的算法。
运算服务器102还可以存在中值滤波法，该方法可以对已经成像的地质雷达信号进行去噪，例如在地质雷达信号传输过程中收到椒盐噪声、高斯噪声的影响，为了消除这两种噪声，改善地质雷达信号的图像，通过增强图像的质量，对该图像直接进行去噪，通过对图像内每一个像素点的灰度值设置为该点某邻域窗口内的所有象素点灰度值的中值，进而平滑图像，达到去噪的效果，该邻域窗口的形状可以为线状、方形、圆形和十字形等，但是中值滤波仅仅只能对地质雷达信号对应的图像直接进行处理，所以计算过程相对复杂，且消除噪声的精度并不高。
运算服务器102还可以存在F-K滤波法，F-K滤波法在地质探勘中是为了消除信号中的噪声，并最大限度的保留有效反射波，可以根据各类噪声和有效反射波的特点进行处理，一般来说在频域中进行滤波去噪时，无法完美的压制有效波，增强干扰波，所以在实际生产中主要以F-K滤波为主的二维的频率波数域为主的二维滤波来压制干扰波，通过二维傅里叶变换将地震数据从T-X变换至F-K，利用有效波和干扰波在传播方向上的差异性，从而使有效波和干扰波在频率波数域分开，进而保证在不伤害有效波的情况下将干扰波剔除，因为F-K滤波法需要在二维傅里叶变换的基础上进行，所以F-K滤波法存在于傅里叶变换算法相同的缺陷，即需要提取所有信号，无法进行局部去噪。
运算服务器102还可以存在本文的地质雷达信号去噪方法，可以将局部的地质雷达信号或者全频段的地质雷达信号进行去噪，并在去噪声过程中，接收到来自控制终端103选定的斜率比值，或者运算服务器102可以进行自动模式，自动的为用户选定斜率比值最大对应的特征值，然后根据该斜率比值筛选噪声和有效信号，将噪声去除后，将有线信号重新构建并得到用户所需要的去噪后的地质雷达信号。
需要说明的是，在运算服务器102中存在的何种去噪算法，可以根据用户通过控制终端103进行切换，还可以默认使用本文的地质雷达信号去噪方法，本文对于此不做限定。
需要说明的是，本文中的运算服务器102可以是与控制终端103一体的，也可以是分立的，运算服务器102可以与地质雷达101有线连接，也可以通过无线网络进行连接，本文对地质雷达101的连接方式不做限定，且该地质雷达101不仅可以进行地质雷达信号探测，也可以进行水声雷达信号探测或者空气介质雷达信号探测，本文对于地质雷达101的属性不做限定。
在地质雷达去噪过程中，受到现有技术中计算复杂、参数冗余、适用性差且无法有效的抑制地质雷达信号噪声的问题，所以亟需一种方法，可以快速准确地去除地质雷达信号中的噪声。
为了解决上述问题，本文实施例提供了一种地质雷达信号去噪方法，能够建立预测砂体厚度的拟合方程，图2是本文实施例提供的一种地质雷达信号去噪方法的步骤示意图，可以应用于信号处理技术领域，可用于地质勘探领域，使用该方法可以解决上述问题，本说明书提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的劳动可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的系统或装置产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行。具体的如图2所示，所述方法可以包括：
步骤201、根据地质雷达信号，确定原始矩阵。
步骤202、计算所述原始矩阵的特征值、特征向量及特征值间的斜率比值。
步骤203、根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵。
步骤204、根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号。
需要说明的是，地质雷达信号存在少量的噪声，但是该噪声的总值远远小于有效信号，将地质雷达信号确定为原始矩阵是较为常规的技术手段，在现有数学工具中，可以将地质雷达信号使用MATLAB-SPAMS稀疏建模工具箱，将地质雷达信号进行稀疏建模并进行矩阵化，该矩阵化的方法存在很多种，本领域技术人员可以根据实际需要更改矩阵化方法，本文在此不做限定。
而特征值和特征向量可以通过特征多项式得到，具体计算过程参见后续实施例，此处不再详述。
在本文中去噪矩阵对应有效信号，即用户期望获得的没有噪声影响的信号所组成的矩阵。
通过上述方法，可以实现对原始地质雷达信号进行出噪，且去噪过程需要的参数简单，计算过程简单，并可以对任何频段的地质雷达进行去噪，还可以局部去噪，无需获取全部的地质雷达信号。
作为本文的一个实施例，步骤202所述计算所述原始矩阵的特征值和特征向量，进一步包括：
对所述原始矩阵进行零均值化处理，得到零均值化矩阵。
根据所述零均值化矩阵，确定所述协方差矩阵。
根据所述协方差矩阵，计算所述原始矩阵的特征值和特征向量。
在本步骤中，零均值化处理可以降低矩阵运算的难度，并且可以加快矩阵在一些运算的收敛速度，例如在数组[3，4，5，6]中，将该数组进行零均值化处理后可以得到[-1.5，-0.5，0.5，1.5]，通过这种方式可以大大加快运算速度，实现本文对地质雷达信号去噪的快速响应，
根据公式Ynm＝Xnm-DXnm(1/n)，可以计算得到中间变量Ynm，其中Xnm为零均值化矩阵，D为单位矩阵，将Xnm以矩阵表示为
Figure BDA0003364280970000081
该零均值化矩阵为n行m列组成，这样可以将公式Ynm＝Xnm-EXnm(1/n)矩阵化表示为
Figure BDA0003364280970000082
单位矩阵E为
Figure BDA0003364280970000091
将中间变量Ynm进行转置，得到中间变量Ynm的转置矩阵
Figure BDA0003364280970000092
对于本领域的技术人员来说将一个矩阵转置得到转置矩阵是很容易实现的，所以本文在此不再赘述
通过公式
Figure BDA0003364280970000093
即可得到协方差矩阵Dmm，为了方便说明，可以将协方差矩阵Dmm使用矩阵表示为
Figure BDA0003364280970000094
使用特征多项式|Dmm-λEmm|计算特征值，然后将特征值带入到齐次线性方程组，即可得到每一个特征值对应的特征向量。
作为本文的一个实施例，步骤所述对所述原始矩阵进行零均值化处理，得到零均值化矩阵，进一步包括：
根据所述原始矩阵中每一列的平均值，构建平均值矩阵。
利用所述原始矩阵减去所述平均值矩阵，得到所述零均值化矩阵。
在本步骤中，计算每一列的平均值，构建平均值矩阵是为了降低每一列，或者每一维的差距，将每一列的平均值计算并得到
Figure BDA0003364280970000095
矩阵化形式为
Figure BDA0003364280970000096
在本文中地质雷达信号的原始矩阵为dnm，dnm的矩阵表达形式为
Figure BDA0003364280970000097
根据上述的公式得到
Figure BDA0003364280970000098
其中Xnm为零均值化矩阵，为了方便说明，可以将零均值化矩阵表示为
Figure BDA0003364280970000101
作为本文的一个实施例，步骤所述根据所述协方差矩阵，计算所述原始矩阵的特征值和特征向量，进一步包括：
根据所述协方差矩阵、单位矩阵和零向量确定所述原始矩阵的特征值；
将所述特征值依次带入至齐次线性方程组中，计算得到每个所述特征值对应的特征向量，其中，所述齐次方程组由所述协方差矩阵、所述单位矩阵和所述零向量组成。
在本步骤中，特征多项式|Dmm-λEmm|中λ为特征值，Emm为单位向量，将|Dmm-λEmm|矩阵化为
Figure BDA0003364280970000102
求解特征多项式|Dmm-λEmm|＝0的所有根λi，其中i＝1，2，3……，在本公式中0表示零向量，并未数字0，λi表示第i个特征值。
根据协方差矩阵Dmm、单位矩阵Emm和零向量0，构建齐次线性方程组(Dmm-λEmm)x＝0，将所有的特征值λ带入到齐次线性方程组中，可以得到每个特征值λi对应的特征向量Vi，将特征向量Vi以矩阵形式表示为
Figure BDA0003364280970000103
作为本文的一个实施例，在步骤203所述计算特征值间的斜率比值，包括：
将所述特征值按照由大到小的顺序进行排序。
将所述特征值的排序序号作为横坐标，将特征值作为纵坐标，建立坐标系。
在所述坐标系下，计算相邻坐标点之间的斜率；
根据所述相邻特征值间的斜率，计算相邻斜率间的斜率比值，将相邻斜率间的斜率比值作为排序后特征值间的斜率比值。
需要说明的是，地质雷达信号中包括有效信号和噪声，噪声是干扰用户获取地质信息的信号，而有效信号是用户期望得到的信号，特征值的大小代表了该特征值对应的特征向量占得地质雷达信号中的权重，特征值越大，在地质雷达信号中起的作用越大。理论上，噪声的占比非常小，所以将所有的特征值进行排序，一般越小的特征值，其特征向量表示噪音。
需要说明的是，计算斜率的目的是为了找出噪声，因为代表噪声的特征向量很小，从有效信号的特征向量过渡到噪声时，会出现一个突变，由大值到小值，变化较大，即斜率会很大，通过找到这个大斜率点，即斜率比值最大的突变点，就能找到有效信号和噪音的分割位置。
在本步骤中，将特征值进行排序，例如数组中[λ1＝3，λ2＝1，λ3＝5，λ4＝6，λ5＝7]，那么按照从大到小的顺序，将数组排序为A＝[λ5，λ4，λ3，λ1，λ2]，然后重新进行排序，确定数组Qi＝[λ1，λ2，λ3，λ4，λ5]，并进行重新赋值将A的至赋给Qi，此时Qi中[λ1＝7，λ2＝6，λ3＝5，λ4＝3，λ5＝1]，然后以当前的Qi中λi的角标作为横坐标，将对应的特征值作为纵坐标，画出坐标系后，根据相邻坐标点间的斜率，例如λ1和λ2之间的斜率为1，λ2和λ3之间的斜率为1，λ3和λ4之间的斜率为2，λ4和λ5之间的斜率为3，然后根据公式
Figure BDA0003364280970000111
计算相邻斜率之间的比值，例如在λ3和λ4，与λ4和λ5之间的斜率进行比较，那么将该斜率比值作为公式中λi的斜率比值即λ5，也就是三个点λ5，λ4，λ3间，角标最大的点的斜率比值。
作为本文的一个实施例，步骤203所述根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵，进一步包括：
根据所述原始矩阵的特征值之间的比值关系确定分界点，并根据所述分界点区分所述原始矩阵中的噪声和有效信号。
将所述有效信号对应的特征向量组合确定所述去噪矩阵。
步骤根据所述原始矩阵的特征值之间的比值关系确定分界点，并根据所述分界点区分所述原始矩阵中的噪声和有效信号，具体包括：
选取最大斜率比值相关的特征值中最大的排序序号。
将特征值排序中位于选取出的排序序号之前的特征值作为所述去噪特征值。
将特征值排序中位于选取出的排序序号之前的特征值作为所述去噪特征值。具体为根据所述原始矩阵的特征值之间的比值关系确定分界点，并根据所述分界点区分所述原始矩阵中的噪声和有效信号，具体为在地质雷达信号中，噪声是相对于有效信号较少的，所以其对于整个地质雷达信号的贡献率是相对较少的，在上述坐标系中，根据斜率的比值，可以表征出有效信号和噪声对于地质雷达信号的贡献率，根据贡献率划定分界点，当在该坐标系中，斜率比值达到最大值时，表征该点前皆为有效信号，而之后皆为噪声。
在本步骤中，选取所述原始矩阵中斜率比值最大的特征值的排序序号。
根据上述Qi＝[λ1，λ2，λ3，λ4，λ5]可以计算出，当选择λ5作为分界点时，该斜率比值是最大的，也就是说，λ1，λ2，λ3，λ4，λ5皆是有效信号，当然，在实际情况中不会出现这样的现象，当然该例子只作为参考，本领域技术人员可以根据实际情况，利用斜率筛选出有效信号和噪声，且在地质雷达信号中一定存在噪声。
步骤将所述有效信号对应的特征向量组合确定所述去噪矩阵，具体包括：
将所述去噪特征值对应的特征向量组合确定所述去噪矩阵。
将选取出来的特征值λi对应的
Figure BDA0003364280970000121
进行重新组合，得到去噪矩阵PnK，将PnK使用矩阵化表示为
Figure BDA0003364280970000122
其中PnK中特征值的顺序已经被按照特征值的大小顺序被打乱，所以需要重新排列，其中PnK的角标K对应选出的特征值的角标，例如本文例子中特征值角标为5，那么K＝5。
作为本文的一个实施例，步骤204所述根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号，进一步包括：
根据所述原始矩阵中每一列的平均值，构建平均值矩阵。
利用所述原始矩阵减去所述平均值矩阵，得到所述零均值化矩阵。
将所述去噪矩阵和所述零均值化矩阵相乘，得到降维矩阵。
对所述去噪矩阵进行转置处理，得到转置去噪矩阵。
将所述转置去噪矩阵和所述降维矩阵相乘，再与所述平均值矩阵相加，确定所述去噪地质雷达信号。
因为去噪矩阵PnK已经被降维到K阶，所以利用PnK与归零化矩阵Xnm相乘，可以得到中间变量YKm，因为PnK是不确定的，所以本文不演示PnK与Xnm相乘的过程，本领域技术人员可以根据实际情况，进行矩阵相乘。
本文的中间变量YKm的矩阵形式为
Figure BDA0003364280970000131
然后利用公式
Figure BDA0003364280970000132
进行重构，需要说明的是
Figure BDA0003364280970000133
为PnK的转置矩阵，本领域技术人员可以很容易得到，而
Figure BDA0003364280970000134
前文已经得到，所以
Figure BDA0003364280970000135
去噪地质雷达信号是很容易得到，
Figure BDA0003364280970000136
的矩阵形式为
Figure BDA0003364280970000137
通过上述方式，可以实现对任意频段的地质雷达信号进行去噪，且去噪过程简单，所需参数简洁，并且在进行分界时，还可以根据用户需要，改变分界条件，更加精细的进行去噪。
如图3所示一种地质雷达信号去噪装置的示意图，包括：
矩阵化单元301，用于根据所述地质雷达信号，确定原始矩阵。
分析单元302，用于计算所述原始矩阵的特征值、特征向量及特征值间的斜率比值。
筛选单元303，用于根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵。
重构单元304，用于根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号。
通过上述装置，可以快速的确定地质雷达信号对应的矩阵，并提取矩阵的特征值和特征向量，并可以根据特征值之间的斜率比值，区分出噪声和有效信号，将有效信号重构，得到去噪的地质雷达信号，且该装置获取到的矩阵可以是全频段的，也可以是局部的信号，并且本装置对于噪声并无要求，对所有噪声皆可进行去噪，普适性强，响应优良。
如图4所示一种地质雷达信号去噪装置的重构单元示意图，作为本文的一个实施例，重构单元304还包括：
平均值模块3041，用于根据所述原始矩阵中每一列的平均值，构建平均值矩阵。
零均值化模块3042，用于利用所述原始矩阵减去所述平均值矩阵，得到所述零均值化矩阵。
降维模块3043，用于将所述去噪矩阵和所述零均值化矩阵相乘，得到降维矩阵。
转置模块3044，用于对所述去噪矩阵进行转置处理，得到转置去噪矩阵。
运算模块3045，用于将所述转置去噪矩阵和所述降维矩阵相乘，再与所述平均值矩阵相加，确定所述去噪地质雷达信号。
如图5所示一种地质雷达信号去噪方法的数据流程图，包括：
步骤501、地质雷达101接收地质雷达信号，并发送至运算服务器102。
步骤502、用户通过控制终端103选定算法，例如本文的地质雷达信号去噪方法，并将指令发送至运算服务器102。
步骤503、运算服务器102对地质雷达信号进行矩阵化，得到原始矩阵。
步骤504、运算服务器102计算原始矩阵的特征值和特征向量。
步骤505、如图6所示一种地质雷达信号去噪方法的相邻特征值间的斜率比值示意图，运算服务器102计算特征值间的斜率比值，并将该斜率比值发送至控制终端103。
步骤506、如图7所示一种地质雷达信号去噪方法的筛选条件示意图，用户可以根据需要选择不同的K值，即筛选条件，也可以输入自动配置，令运算服务器102选择最优解，并进行自动计算，在图中，原始单道记录为地质雷达信号中的一道信号，从不同的PCA对应的K值的拟合曲线可以看出，在本文中K＝5，与原始单道记录即地质雷达信号最为拟合，且将井盖响应即噪声去除。
步骤507、运算服务器102根据K值计算去噪矩阵。
步骤508、运算服务器102根据去噪矩阵和原始矩阵，得到去噪地质雷达信号并发送至控制终端。
通过上述的方式，可以实现切换不同去噪算法，以解决不同类型地质雷达信号去噪问题，还可以通过人机交互的过程，令用户确定K值，进行更加精细化的筛选，而且运算服务器102还具有自动配置功能，可以在计算完成特征值间斜率后，直接得出最优解，并在得到最优解后，向控制终端103发送去噪后的地质雷达信号。
如图8所示，为本文实施例提供的一种计算机设备，所述计算机设备802可以包括一个或多个处理器804，诸如一个或多个中央处理单元(CPU)，每个处理单元可以实现一个或多个硬件线程。计算机设备802还可以包括任何存储器806，其用于存储诸如代码、设置、数据等之类的任何种类的信息。非限制性的，比如，存储器806可以包括以下任一项或多种组合：任何类型的RAM，任何类型的ROM，闪存设备，硬盘，光盘等。更一般地，任何存储器都可以使用任何技术来存储信息。进一步地，任何存储器可以提供信息的易失性或非易失性保留。进一步地，任何存储器可以表示计算机设备802的固定或可移除部件。在一种情况下，当处理器804执行被存储在任何存储器或存储器的组合中的相关联的指令时，计算机设备802可以执行相关联指令的任一操作。计算机设备802还包括用于与任何存储器交互的一个或多个驱动机构808，诸如硬盘驱动机构、光盘驱动机构等。
计算机设备802还可以包括输入/输出模块810(I/O)，其用于接收各种输入(经由输入设备812)和用于提供各种输出(经由输出设备814))。一个具体输出机构可以包括呈现设备816和相关联的图形用户接口(GUI)818。在其他实施例中，还可以不包括输入/输出模块810(I/O)、输入设备812以及输出设备814，仅作为网络中的一台计算机设备。计算机设备802还可以包括一个或多个网络接口820，其用于经由一个或多个通信链路822与其他设备交换数据。一个或多个通信总线824将上文所描述的部件耦合在一起。
通信链路822可以以任何方式实现，例如，通过局域网、广域网(例如，因特网)、点对点连接等、或其任何组合。通信链路822可以包括由任何协议或协议组合支配的硬连线链路、无线链路、路由器、网关功能、名称服务器等的任何组合。
对应于图2和图5中的方法，本文实施例还提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器运行时执行上述方法的步骤。
本文实施例还提供一种计算机可读指令，其中当处理器执行所述指令时，其中的程序使得处理器执行如图2和图5所示的方法。
应理解，在本文的各种实施例中，上述各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本文实施例的实施过程构成任何限定。
还应理解，在本文实施例中，术语“和/或”仅仅是一种描述关联对象的关联关系，表示可以存在三种关系。例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。
本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本文的范围。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
在本文所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另外，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口、装置或单元的间接耦合或通信连接，也可以是电的，机械的或其它的形式连接。
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本文实施例方案的目的。
另外，在本文各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以是两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本文的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本文各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
本文中应用了具体实施例对本文的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本文的方法及其核心思想；同时，对于本领域的一般技术人员，依据本文的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本文的限制。

Claims (10)

1.一种地质雷达信号去噪方法，其特征在于，包括：
根据地质雷达信号，确定原始矩阵；
计算所述原始矩阵的特征值、特征向量及特征值间的斜率比值；
根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵；
根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号。
2.根据权利要求1所述的地质雷达信号去噪方法，其特征在于，所述计算所述原始矩阵的特征值、特征向量，进一步包括：
对所述原始矩阵进行零均值化处理，得到零均值化矩阵；
根据所述零均值化矩阵，确定协方差矩阵；
根据所述协方差矩阵，计算所述原始矩阵的特征值和特征向量。
3.根据权利要求2所述的地质雷达信号去噪方法，其特征在于，所述对所述原始矩阵进行零均值化处理，得到零均值化矩阵，进一步包括：
根据所述原始矩阵中每一列的平均值，构建平均值矩阵；
利用所述原始矩阵减去所述平均值矩阵，得到所述零均值化矩阵。
4.根据权利要求2所述的地质雷达信号去噪方法，其特征在于，所述根据所述协方差矩阵，计算所述原始矩阵的特征值和特征向量，进一步包括：
根据所述协方差矩阵、单位矩阵和零向量确定所述原始矩阵的特征值；
将所述特征值依次带入至齐次线性方程组中，计算得到每个所述特征值对应的特征向量，其中，所述齐次线性方程组由所述协方差矩阵、所述单位矩阵和所述零向量组成。
5.根据权利要求1所述的地质雷达信号去噪方法，其特征在于，计算特征值间的斜率比值，包括：
将所述特征值按照由大到小的顺序进行排序；
将所述特征值的排序序号作为横坐标，将所述特征值作为纵坐标，建立坐标系；
在所述坐标系下，计算相邻坐标点之间的斜率；
根据所述相邻特征值间的斜率，计算相邻斜率间的斜率比值，将相邻斜率间的斜率比值作为排序后特征值间的斜率比值。
6.根据权利要求5所述的地质雷达信号去噪方法，其特征在于，所述根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵，进一步包括：
选取最大斜率比值相关的特征值中最大的排序序号；
将特征值排序中位于选取出的排序序号之前的特征值作为所述去噪特征值；
将所述去噪特征值对应的特征向量组合，确定所述去噪矩阵。
7.根据权利要求1所述的地质雷达信号去噪方法，其特征在于，所述根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号，进一步包括：
根据所述原始矩阵中每一列的平均值，构建平均值矩阵；
利用所述原始矩阵减去所述平均值矩阵，得到零均值化矩阵；
将所述去噪矩阵和所述零均值化矩阵相乘，得到降维矩阵；
对所述去噪矩阵进行转置处理，得到转置去噪矩阵；
将所述转置去噪矩阵和所述降维矩阵相乘，再与所述平均值矩阵相加，确定所述去噪地质雷达信号。
8.一种地质雷达信号去噪装置，其特征在于，包括：
矩阵化单元，用于根据地质雷达信号，确定原始矩阵；
分析单元，用于计算所述原始矩阵的特征值、特征向量及特征值间的斜率比值；
筛选单元，用于根据所述特征值间的斜率比值，从所述原始矩阵的特征值中筛选出去噪特征值，并根据所述去噪特征值确定去噪矩阵；
重构单元，用于根据所述去噪矩阵和所述原始矩阵，确定去噪地质雷达信号。
9.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1-7任意一项所述地质雷达信号去噪方法。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1-7任意一项所述地质雷达信号去噪方法。
CN202111400693.XA 2021-11-19 2021-11-19 一种地质雷达信号去噪方法、装置、设备及介质 Pending CN114154537A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111400693.XA CN114154537A (zh) 2021-11-19 2021-11-19 一种地质雷达信号去噪方法、装置、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111400693.XA CN114154537A (zh) 2021-11-19 2021-11-19 一种地质雷达信号去噪方法、装置、设备及介质

Publications (1)

Publication Number Publication Date
CN114154537A true CN114154537A (zh) 2022-03-08

Family

ID=80457250

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111400693.XA Pending CN114154537A (zh) 2021-11-19 2021-11-19 一种地质雷达信号去噪方法、装置、设备及介质

Country Status (1)

Country Link
CN (1) CN114154537A (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109085547A (zh) * 2018-09-26 2018-12-25 湖南时变通讯科技有限公司 一种表层穿透雷达回波信号的去噪方法和相关装置
CN112472079A (zh) * 2020-11-23 2021-03-12 青岛歌尔智能传感器有限公司 血氧饱和度检测装置、设备及存储介质

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109085547A (zh) * 2018-09-26 2018-12-25 湖南时变通讯科技有限公司 一种表层穿透雷达回波信号的去噪方法和相关装置
CN112472079A (zh) * 2020-11-23 2021-03-12 青岛歌尔智能传感器有限公司 血氧饱和度检测装置、设备及存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
ACCEPTEDLIN: "特征提取之PCA", pages 1 - 13, Retrieved from the Internet <URL:https://blog.csdn.net/u013185349/article/details/78084576> *
沈宏益著: "《生态安全与西藏新型工业化研究》", 30 November 2014, 厦门：厦门大学出版社, pages: 172 - 182 *

Similar Documents

Publication Publication Date Title
Huang et al. Signal extraction using randomized-order multichannel singular spectrum analysis
JP2999266B2 (ja) １次元又は２次元の方向推定または周波数推定のための信号の高分解能評価方法
Qiao et al. Gridless line spectrum estimation and low-rank Toeplitz matrix compression using structured samplers: A regularization-free approach
CN110276726B (zh) 一种基于多通道网络先验信息引导的图像去模糊方法
Boashash et al. Robust multisensor time–frequency signal processing: A tutorial review with illustrations of performance enhancement in selected application areas
EP2819025A1 (en) Method for reducing noise in data-sets of harmonic signals
CN108710150A (zh) 一种基于稳健奇异谱分析的地震不规则噪声去除方法
EP3994690A1 (en) Audio processing apparatus and method for denoising a multi-channel audio signal
CN109506135A (zh) 管道泄漏点定位方法及装置
CN111863014A (zh) 一种音频处理方法、装置、电子设备和可读存储介质
Zhou et al. A hybrid method for noise suppression using variational mode decomposition and singular spectrum analysis
CN115409735A (zh) 一种InSAR干涉图像的滤波方法及装置
CN117111000A (zh) 一种基于双通道注意力残差网络的sar梳状谱干扰抑制方法
CN116699526A (zh) 一种基于稀疏与低秩模型的车载毫米波雷达干扰抑制方法
CN112285793B (zh) 一种大地电磁去噪方法及系统
CN107481732B (zh) 一种口语测评中的降噪方法、装置及终端设备
Li et al. Desert seismic signal denoising by 2D compact variational mode decomposition
CN106022358A (zh) 一种高光谱图像分类方法和装置
US11482239B2 (en) Joint source localization and separation method for acoustic sources
CN110929759B (zh) 检测模型的训练装置、方法及心电数据处理方法、装置
CN114154537A (zh) 一种地质雷达信号去噪方法、装置、设备及介质
CN116612735A (zh) 一种水声音频去噪方法
KR100633555B1 (ko) 신호 복구 방법
Gantenapalli et al. Selective mean filtering for reducing impulse noise in digital color images
CN114047481A (zh) 一种基于子空间正交性的稳健自适应波束形成方法

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