CN112765546B - 一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质 - Google Patents
一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质 Download PDFInfo
- Publication number
- CN112765546B CN112765546B CN202110038837.5A CN202110038837A CN112765546B CN 112765546 B CN112765546 B CN 112765546B CN 202110038837 A CN202110038837 A CN 202110038837A CN 112765546 B CN112765546 B CN 112765546B
- Authority
- CN
- China
- Prior art keywords
- signal
- equivalent
- transmitted
- power frequency
- frequency noise
- 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 84
- 230000009467 reduction Effects 0.000 title claims abstract description 27
- 230000005284 excitation Effects 0.000 claims abstract description 42
- 230000004044 response Effects 0.000 claims description 58
- 230000005540 biological transmission Effects 0.000 claims description 37
- 238000012545 processing Methods 0.000 claims description 12
- 238000005070 sampling Methods 0.000 claims description 9
- 238000005311 autocorrelation function Methods 0.000 claims description 7
- 238000001914 filtration Methods 0.000 claims description 6
- 238000005314 correlation function Methods 0.000 claims description 5
- 238000004891 communication Methods 0.000 claims description 4
- 238000004590 computer program Methods 0.000 claims description 4
- 230000000694 effects Effects 0.000 description 11
- 230000008569 process Effects 0.000 description 9
- 238000011410 subtraction method Methods 0.000 description 5
- 238000013461 design Methods 0.000 description 4
- 230000006870 function Effects 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000008054 signal transmission Effects 0.000 description 1
- 230000009131 signaling function Effects 0.000 description 1
Classifications
-
- 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/15—Correlation function computation including computation of convolution operations
- G06F17/153—Multidimensional correlation or convolution
-
- 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)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Noise Elimination (AREA)
Abstract
本发明公开了一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质，所述方法包括如下步骤：设定发射信号的时长T与工频噪声周期Tc的关系；连续发出发射信号传送至大地系统，并基于T与Tc的关系确定等效发射信号，其中，基于等效发射信号得到的等效接收信号中不包含工频噪声；若T为Tc的整数倍，等效发射信号为连续发射的相邻两组发射信号相减得到，等效接收信号为相邻两组接收信号相减得到；若T为Tc的整数倍加二分之一个Tc，等效发射信号为连续发射的相邻两组发射信号相加得到，等效接收信号为相邻两组接收信号相加得到。本发明针对工频干扰的特点，通过优化激励信号来抵消工频噪声，可以高效去除工频噪声。
Description
技术领域
本发明属于大地电磁技术领域，具体涉及一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质。
背景技术
由于接收端测量的信号极其微弱，可控源电磁法非常容易受到各种噪声的影响。随着工业化的进展，大量的高压输电线和各种用电设施所带来的工频干扰无处不在，目前几乎不存在无电力噪声影响的电磁区域。在工业用电密集区域，工频干扰极大的影响了测量数据的精确性，导致信噪比大大降低，影响了大地系统的辨识精度，最终无法正确的提取大地信息。因此，去除工频干扰是进行地球物理高精度勘探的基本前提。
去除工频干扰传统的方法主要包括：陷波法和工频噪声估计法。陷波法直接在频率域滤除特定频点的干扰，其原理简单，但存在两个较明显的问题：一是滤波器在滤除特定频点干扰的同时也将同一频点的有用信号滤掉了；二是实际的陷波器会导致特定频点附近频率的响应失真。工频噪声估计法，也称工频噪声减去法。减去法包括块减去法和正弦减去法，原理都是将估计出的工频噪声从测量信号中减掉。只不过块减法是通过增加一个单独的噪声记录仪来记录噪声，而正弦减去法则是使用记录的噪声来估计工频干扰的幅值和相位。工频噪声减去法对消除工频干扰有较好的效果，但前提是噪声在短时间内是稳定的，且野外工作时需要增加噪声记录仪，反而增加了野外的工作量，相对较为麻烦且对估计算法要求较高。因此，如何更加便捷以及有效降低噪声干扰，是亟待研究的。
发明内容
本发明的目的是提出一种抵消工频噪声的技术手段，进而提出了一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质，所述方法不同于现有去除工频噪声的技术手段，其针对工频干扰的特点，通过优化激励信号来抵消工频噪声，可以高效去除工频噪声，为后续辨识大地系统辨识奠定了基础。
一方面，本发明提供一种应用于可控源电磁法相关辨识中的降噪方法，包括如下步骤：
设定发射信号的时长T与工频噪声周期Tc的关系；
连续发出发射信号传送至大地系统，并基于T与Tc的关系确定等效发射信号，其中，基于等效发射信号得到的等效接收信号中不包含工频噪声；
若发射信号的时长T为工频噪声周期Tc的整数倍，等效发射信号为连续发射的相邻两组发射信号相减得到，等效接收信号为相邻两组接收信号相减得到；若发射信号的时长T为工频噪声周期Tc的整数倍加二分之一个Tc，等效发射信号为连续发射的相邻两组发射信号相加得到，等效接收信号为相邻两组接收信号相加得到。
本发明基于工频噪声在短时间内是稳定的前提下，考虑到工频噪声具有周期性和对称性的特点，总结出上述两种激励信号优化方式来抵消工频噪声，为后续辨识大地系统辨识奠定了基础。
可选地，所述发射信号的时长T等于发射信号的有效激励时间与零激励时间之和。
本发明消除工频噪声的前提是连续发射两次激励信号，且两激励信号的响应相互不影响，即保证每次激励信号输入后要充分的获得各自的响应。故而要求每次在有效激励信号发射之后要有一段零激励时间，这段零激励时间不仅要保证发射信号充分的激励大地系统，使接收端获得完整的激励响应，同时，零激励时间加上有效激励时间作为发射信号的时长T，其应该满足与工频噪声周期的特定关系。因此，零激励时间等于发射信号的时长T减去发射信号的有效激励时间。
可选地，若发射信号为m序列信号，发射信号的有效激励时间以及零激励时间对应的补零的个数如下所示：
Tv＝(Np*Lm*Nspb)/fs
式中，Npz为零激励时间对应的补零的个数，Tpz为零激励时间，Tv为发射信号的有效激励时间，Np为m序列的周期数，Lm为一个m序列周期的码元数，Nspb为每个码元的采样数，fs为采样频率，K为正整数；
其中，case1表示发射信号的时长T为工频噪声周期Tc的K倍，case2表示发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc。
可选地，还包括利用等效发射信号和等效接收信号辨识大地脉冲响应，其中，等效发射信号需满足：
case1：当发射信号的时长T为工频噪声周期Tc的K倍时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝-xm(t)，对应得到的等效发射信号为m序列2xm(t)；
case2：当发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝xm(t)，对应得到的等效发射信号为m序列2xm(t)，其中，xm(t)为m序列。
在实际情况中，工频噪声几乎不可能是理想的信号，它的幅值和相位可能有轻微的畸变。这种情况下的工频噪声可以视为理想工频噪声加上随机噪声。为了消除随机噪声，设置m序列，利用相关辨识方法，在求解大地脉冲响应的过程中，可将随机噪声消除掉。
可选地，所述辨识大地脉冲响应是利用m序列自相关特点采用维纳滤波反卷积法计算得到，其中，大地脉冲响应ht如下所示：
其中，为等效接收信号与等效发射信号的互相关函数，/>为等效发射信号的自相关函数；采用反卷积求解大地脉冲响应ht存在：
其中，Aτ是等效发射信号st的归一化自相关值，Bτ是等效发射信号st和等效接收信号Vt的归一化互相关值，τ＝0,1,…,n；n为采样点数，根据奈奎斯特-香农采样定理可知，当时间信号函数ht的最高频率为fH时，采样频率必须满足fs≥2fH,一般取(3～5)·fH，故采样点数一般为n＝(3～5)·fH，具体取值应根据计算量和求解精度来确定。Aτ、Bτ分别表示如下：
其中，是等效接收信号Vt的自相关函数在零延时处的相关值，且等效接收信号Vt的自相关函数/>、等效接收信号与等效发射信号的互相关函数/>、等效发射信号的自相关函数/>均利用现有技术可以计算得出。等效发射信号st和等效接收信号Vt分别表示如下：
y1(t)和y2(t)为连续发射的相邻两组发射信号x1(t)、x2(t)对应的响应信号。
可选地，case1：当发射信号的时长T为工频噪声周期Tc的K倍时，等效发射信号的功率为设计的相邻两段发射功率总和；case2：当发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc时，等效发射信号的功率等于设计的发射功率。
相邻两组发射信号x1(t)、x2(t)在相加的过程中，其等效响应得到加强，其功率是原两段系统响应功率之和，进而等效发射信号也得到加强，其等效发射功率是原发射功率的总和。这使得后续若进行系统辨识求解时计算误差可以减小，提高了系统辨识求解精度。
可选地，若信号的发射端或响应信号的接收端存在固有误差，选择发射信号的时长T为工频噪声周期Tc的整数倍，等效发射信号为连续发射的相邻两组发射信号相减得到。
如果发射端或接受端存在某种固有的误差，比如说仪器固有误差或者外界直流干扰等情况，相减抵消法则可以在抵消工频干扰的情况下，同时抵消掉这类固有误差。若采用相加的方法，反而会将固有误差叠加放大，导致辨识精度相对低一些。
第二方面，本发明提供一种应用于可控源电磁法相关辨识中的降噪系统，包括发送端以及接收端，其中，发送端连续发出发射信号传送至大地系统，接收端接收经过大地系统传输后的响应信号；
若发射信号的时长T为工频噪声周期Tc的整数倍，连续发射的相邻两组发射信号相减得到等效发射信号；若发射信号的时长T为工频噪声周期Tc的整数倍加二分之一个Tc，连续发射的相邻两组发射信号相加得到等效发射信号，其中，基于等效发射信号对应的等效接收信号中不包含工频噪声。
可选地，所述降噪系统，还包括处理模块，用于利用等效发射信号和等效接收信号辨识大地脉冲响应，其中，等效发射信号需满足：
case1：当发射信号的时长T为工频噪声周期Tc的K倍时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝-xm(t)，对应得到的等效发射信号为m序列2xm(t)；
case2：当发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝xm(t)，对应得到的等效发射信号为m序列2xm(t)，其中，xm(t)为m序列；
所述处理模块设置于所述接收端内，或与所述接收端外并与所述接收端通信连接。
其中，处理模块可以是内置于接收端的处理器等具有数据处理功能的元件，也可以是与接收端通信连接的外设电子终端，譬如，计算机等。
第三方面，本发明提供一种可读存在介质，存储了计算机程序，所述计算机程序被处理器调用以执行所述一种应用于可控源电磁法相关辨识中的降噪方法的步骤。
有益效果
本发明为了抵消工频噪声，利用工频噪声具有周期性的特点，设置发射信号的时长T为工频噪声周期Tc的整数倍，进而将连续发射的相邻两组发射信号相减构建等效发射信号，等效发射信号对应的等效接收信号中工频噪声相抵消，实现了工频去噪；同理，利用工频噪声具有周期性和对称性的特点，设置发射信号的时长T为工频噪声周期Tc的整数倍加二分之一个Tc，进而将连续发射的相邻两组发射信号相加构建等效发射信号，等效发射信号对应的等效接收信号中工频噪声相抵消，同样实现了工频去噪。从而，本发明通过优化激励信号实现去除工频噪声，为后续辨识大地系统辨识奠定了基础。相较于现有的陷波法和工频噪声估计法，更具优势。
在本发明的进一步优选方案中，还考虑到了随机噪声，进而从等效信号方面设计激励信号，设计等效发射信号为m序列，利用m序列的自相关性消除随机噪声，进而辨识出的大地脉冲响应中消除了工频噪声和随机噪声。
附图说明
图1是大地系统辨识简化模型；
图2是T为Tc整数倍的关系图，其中，虚线部分表示有多个整数个周期的正弦信号；
图3是T为Tc的非整数倍的关系图，其中，虚线部分表示有多个整数个周期的正弦信号；
图4是case1情况下的发射信号与其等效信号；
图5是case2情况下的发射信号与其等效信号；
图6是大地系统辨识简化另一模型；
图7是case1情况下的实际发射信号与其等效信号；
图8是case2情况下的实际发射信号与其等效信号；
图9为case1均匀半空间辨识结果与辨识误差，其中，(a)图为辨识结果比对图，(b)图为辨识结果误差图；
图10为case2均匀半空间辨识结果与辨识误差，其中，(a)图为辨识结果比对图，(b)图为辨识结果误差图。
具体实施方式
本发明提供一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质。本发明提供的所述降噪方法至少可以抵消工频噪声，且优选还抵消掉随机噪声，下述将对其原理进行说明。
大地系统辨识模型可以简化成图1。x(t)为发射端的电流信号，h(t)是大地系统的脉冲响应，n(t)为加性噪声，y(t)是接收端的测量电压信号。大地系统可以用以下数学关系式表示：
y(t)＝x(t)*h(t)+n(t) (1)
其中，符号“*”表示卷积。从上述公式(1)可知，接收信号y(t)是由有用信息和噪声组成的。降噪的目的就是消除接收的电压信号中的噪声，最终为辨识出大地脉冲响应奠定基础，提高辨识结果的准确性。而辨识结果的好坏跟接收端的信噪比有直接关系，噪声去除的越干净，辨识的精度就越高。本发明要考虑强工频干扰和随机噪声干扰的消除问题。
关于强工频干扰。
发射信号x(t)由两组序列组成的，分别用x1(t)和x2(t)表示，假设噪声项只有单一的工频噪声，且先不考虑谐波的影响。则有：
y1(t)＝x1(t)*h(t)+nc(t) 0＜t≤T (2)
y2(t+T)＝x2(t+T)*h(t+T)+nc(t+T) 0＜t≤T (3)
信号x1(t)和x2(t)的发射时间长度均为T，y1(t)和y2(t)分别对应x1(t)和x2(t)的含噪系统响应。若要高精度地求解大地系统的单位脉冲响应，根据式(2)和式(3)可知，能否有效去除噪声项是问题的关键所在。下面从每组发射信号的时长T与工频噪声周期的关系方面进行考虑。
发射信号时长T与工频噪声周期Tc之间的关系考虑以下两种情况：
case1:T是Tc的整数倍。
如图2所示，Tc为工频噪声的一个周期，T是Tc的整数倍，即有
T＝KTc K∈N* (4)
其中，T为一次发射的时长，K为正整数。由于工频噪声具有周期性，只需要保证第一次发射信号和第二次发射信号是无间隙的，即x1(t)和x2(t)是连续发射，则发射起始时间可以是任意时刻，故有：
nc(t+T)＝nc(t+KTc)＝nc(t) t≥0 (5)
则式(2)减去式(3)有：
y1(t)-y2(t+T)＝x1(t)*h(t)-x2(t+T)*h(t+T) (6)
只要保证T大于大地系统单位脉冲响应调整时间，大地系统单位脉冲响应调整时间为大地脉冲响应从时刻0到响应值等于0的时刻的时长，x1(t)和x2(t)即可当做两次独立同时间段的发射信号，y1(t)和y2(t)为对应的接收响应。则(6)式可简化为：
y1(t)-y2(t)＝[x1(t)-x2(t)]*h(t) (7)
本发明将x1(t)-x2(t)作为等效发射信号，y1(t)-y2(t)作为等效接收信号。
case2:T是Tc的非整数倍。
如图3所示，Tc为工频噪声的一个周期，T是Tc的整数倍加二分之一个Tc，即有
其中，T为一次发射的时长，K为正整数。在x1(t)和x2(t)是连续发射的前提下，发射起始时间可以是任意时刻，根据工频噪声的周期性和对称性有：
则式(2)加上式(3)有：
y1(t)+y2(t+T)＝x1(t)*h(t)+x2(t+T)*h(t+T) (10)
在保证T大于大地系统单位脉冲响应调整时间的前提下，x1(t)和x2(t)可当做两次独立同时间段的发射，y1(t)和y2(t)为对应的接收响应。则(10)式可简化为：
y1(t)+y2(t)＝[x1(t)+x2(t)]*h(t) (11)
本发明将x1(t)+x2(t)作为等效发射信号，y1(t)+y2(t)作为等效接收信号。
综上所述，case1和case2两种情况是从激励信号的发射时长方面进行了考虑，以上抵消工频噪声的方法是基于工频噪声在短时间内是稳定的前提下推导的。
实际情况中还要考虑残留噪声对系统辨识的影响，上述部分假设噪声是单一的理想工频噪声，但是在实际情况中，工频噪声几乎不可能是理想的信号，它的幅值和相位可能有轻微的畸变。这种情况下的工频噪声可以视为理想工频噪声加上随机噪声。即有：
nc(t)＝ni(t)+nr(t) (12)
其中nc(t)为畸变的工频噪声，ni(t)为理想的工频噪声，nr(t)为随机噪声。通过本文所提消除工频噪声的方法，ni(t)可以完全被抵消，则噪声只剩下随机噪声。下面从发射信号与等效发射信号的关系方面进行考虑，以消除残留噪声、提高系统辨识精度。
在式(7)和式(11)中的等效发射信号[x1(t)-x2(t)]和[x1(t)+x2(t)]，设计其为m序列，这样通过相关辨识求解大地脉冲响应的过程中，就可将随机噪声消除掉。故发射信号的设计原则就是使等效发射信号为m序列，且发射信号本身易于实现。由于该部分是现有技术可实现过程，本发明对其不进行具体的阐述。
下面针对等效发射信号的两种情况，根据发射信号设计原则分别给出两种简单实用的发射信号，具体实现步骤如下。
(a)针对等效发射信号为[x1(t)-x2(t)]的情况，如下：
确定需要的m序列xm(t)。生成m序列的初始状态、反馈系数、阶数、码元宽度、幅值等参数可根据需要自行选定；
令x1(t)＝xm(t)、x2(t)＝-xm(t)。这样等效发射信号[x1(t)-x2(t)]就为m序列2xm(t)。
(b)针对等效发射信号为[x1(t)+x2(t)]的情况，如下：
确定需要的m序列xm(t)。生成m序列的初始状态、反馈系数、阶数、码元宽度、幅值等参数可根据需要自行选定；
令x1(t)＝xm(t)、x2(t)＝xm(t)。这样等效发射信号[x1(t)+x2(t)]就为m序列2xm(t)。
如图4，上部分为设计的发射信号，由xm(t)和-xm(t)两部分组成，下部分为等效发射信号2xm(t)，等效信号的功率为设计发射功率的总和。
如图5，上部分为设计的发射信号，由xm(t)和xm(t)两部分组成，下部分为等效发射信号2xm(t)，等效信号的功率与设计发射功率相等。
以上发射信号的设计仅从等效信号设计成m序列这一原则进行考虑，即仅考虑在抵消工频噪声后，利用m序列的自相关性消除随机噪声。
在进行工频噪声去除处理之后，仍然会有残留的随机噪声。由于等效发射信号为m序列，则可利用其相关特性通过维纳滤波反卷积法对大地脉冲响应进行辨识。假设处理工频噪声之后的系统等效为下图6：
图6中nr为残留的随机噪声，st和Vt分别为等效发射信号和等效接收信号。则各量之间存在如下关系：
Vt＝st*ht+nr (13)
对式(13)两边同时与输入信号st进行相关运算，得到：
(14)
式中，为等效接收信号与等效发射信号的互相关函数，/>为等效发射信号的自相关函数，/>为随机噪声与等效发射信号的互相关函数。利用m序列良好的自相关，等效发射信号是m序列时，其与随机噪声的相关性为零，即/>项几乎为零。式(14)简化为：
(15)
针对式(15)进行反卷积求解大地系统脉冲函数，Anton Ziolkowski(2013)推导并提出了维纳滤波算法：
其中，ht是大地脉冲响应的估计值，Aτ是等效发射信号st的归一化自相关值，Bτ是等效发射信号st和等效接收信号Vt的归一化互相关值。其中，等效发射信号st和等效接收信号Vt分别表示如下：
式(16)中左侧方阵中的元素Aτ和右侧矩阵中的元素Bτ分别为：
通过式(16)～式(20)，可以求出大地脉冲响应，且这种利用m序列自相关特点进行求解的方法，其算法本身具有很强的抗噪性能。这样工频噪声中夹杂着的随机噪声就能够很好的消除。
而抵消工频噪声的过程中，必须要考虑发射周期T与工频噪声周期Tc之间的特定关系，当二者不满足特定关系时，则需要在原有效发射信号之后进行补零，以达到抵消工频噪声的前提条件。具体的补零原则和方法在下文进行详细讨论。
由于本文消除工频噪声的前提是连续发射两次激励信号，且两激励信号的响应相互不影响，即保证每次激励信号输入后要充分的获得各自的响应。这就要求每次在有效激励信号发射之后要有一段零激励时间，这段时间不仅要保证发射信号充分的激励大地系统，使接收端获得完整的激励响应，同时，零激励时间加上有效激励时间应该满足与工频噪声周期的特定关系。
设发射有效信号时间为Tv、补零时间为Tpz，则二者与一次发射总时间T的关系为：
T＝Tv+Tpz (21)
而根据前面推导的发射时间和工频噪声周期之间的关系，式(21)可变为：
其中，case1对应T是Tc的整数倍情况；case2对应T是Tc非整数倍情况。其中Tv的求解如下：
Tv＝(Np*Lm*Nspb)/fs (23)
其中，Np为m序列的周期数，Lm为一个m序列周期的码元数，Nspb为每个码元的采样数，fs为采样频率。
补零的时间可通过式(22)求得，补零的个数Npz可由下式确定：
其中，参数K要保证Npz大于零。
结合发射信号关于等效信号设计部分，补零后实际发射信号及其等效发射信号如附图7和图8所示。
基于上述原理性说明，本发明提供的一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质将结合实施例对其进行进一步的说明。
实施例1：
本实施例提供的一种应用于可控源电磁法相关辨识中的降噪方法，用于抵消工频噪声，其包括如下步骤：
设定发射信号的时长T与工频噪声周期Tc的关系；
连续发出发射信号，并基于T与Tc的关系确定等效发射信号，其中，基于等效发射信号得到的等效接收信号中不包含工频噪声；
case1：若发射信号的时长T为工频噪声周期Tc的整数倍，等效发射信号为连续发射的相邻两组发射信号相减得到。
case2：若发射信号的时长T为工频噪声周期Tc的整数倍加二分之一个Tc，等效发射信号为连续发射的相邻两组发射信号相加得到。
应当理解，本实施例中等效接收信号中不包含工频噪声，实现了工频噪声的消除。为了辨识出大地脉冲响应，还包括基于等效接收信号辨识出大地脉冲响应。其中，基于消除工频噪声的等效接收信号如何辨识出大地脉冲响应的技术手段，本实施例对其不进行具体的限定，任意一种从消除了工频噪声的等效接收信号中辨识出大地脉冲响应的技术手段均属于本实施例的保护范围。
实施例2：
在实施例1的基础上，本实施例还要消除随机噪声，提供的一种应用于可控源电磁法相关辨识中的降噪方法，用于抵消工频噪声，其包括如下步骤：
获取发射信号的时长T与工频噪声周期Tc的关系；
连续发出发射信号，并基于T与Tc的关系确定等效发射信号；
利用等效发射信号和等效接收信号并基于利用m序列自相关特点辨识大地脉冲响应；
case1：若发射信号的时长T为工频噪声周期Tc的整数倍，等效发射信号为连续发射的相邻两组发射信号相减得到；以及连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝-xm(t)，对应得到的等效发射信号为m序列2xm(t)；
case2：若发射信号的时长T为工频噪声周期Tc的整数倍加二分之一个Tc，等效发射信号为连续发射的相邻两组发射信号相加得到；以及连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝xm(t)，对应得到的等效发射信号为m序列2xm(t)，其中，xm(t)为m序列。
本实施例中，按照公式(13)-公式(15)的逻辑可以消除随机噪声，其中，如何根据公式(15)辨识出大地脉冲响应的技术手段，本实施例对其不进行具体的限定。
实施例3：
本实施例在实施例2的基础上，选用反卷积求解大地系统脉冲函数。具体利用上述公式(16)～(20)进行求解公式(15)，计算所得的大地脉冲响应消除了工频噪声和随机噪声。
在实施例1-实施例3的基础上，考虑到零激励时间加上有效激励时间应该满足与工频噪声周期的特定关系。因此，按照公式(21)至公式(24)确定补零的个数Npz。
综上所述，可知，对于case1情形，通过设置每段发射信号的时长为工频噪声周期时长的整数倍，对两段发射信号相减得到等效发射信号，对两段系统响应相减得到等效系统响应(等效接收信号)。在对两段系统响应相减的过程中，工频噪声得到抵消，接收端的信噪比得到了大大的提升。
在工频噪声得到抵消后，进一步将大地系统简化为一个只含有随机噪声的等效系统。通过将等效发射信号设计成m序列，利用m序列的良好自相关性消除随机噪声的影响，最后通过维纳滤波反卷积方法求得大地脉冲响应。去噪效果和辨识效果及辨识误差如附图9所示，从图可知，相减去工频噪声方法去噪效果好，系统辨识精度高。
对于case2情形，通过设置每段发射信号时长是工频噪声周期时长整数陪加二分之一个工频噪声周期的方法，利用工频噪声的周期性和对称性性质，采取对两段系统响应相加的方法抵消掉工频噪声。在相加的过程中，等效响应得到加强，其功率是原两段系统响应功率之和。等效发射信号也得到加强，其等效发射功率是原发射功率的总和。这使得后面相关运算和反卷积运算的计算误差大大减小，极大的提高了系统辨识求解精度。去噪效果和辨识效果及辨识误差如附图10所示，从图可知，相加去工频噪声方法较相减去工频噪声方法去噪效果更好，辨识精度更高。
如果发射端或接受端存在某种固有的误差，比如说仪器固有误差或者外界直流干扰等情况，相减抵消法则可以在抵消工频干扰的情况下，同时抵消掉这类固有误差。对于此类固有误差，case2采用相加的方法，只会将固有误差叠加放大，最终导致辨识精度下降。从消除固有误差方面来讲，case1显然要比case2的去噪效果好。
所以，无论是case1相减去除工频噪声还是case2相加去除工频噪声，都有各自的优缺点，应用场景不同，所得到的效果会有所区别，辨识精度也随之变化。本发明的两种抵消工频噪声的方法对于工频噪声和畸变的工频噪声去噪效果是良好的，去噪方法极大的提高了大地系统的辨识精度，为进一步解释、提取大地信息提供了可靠的依据。本文所提方法原理简单、易于实现、去除特定噪声效果较好，为消除工频噪声提供了一种新的研究思路。
基于上述降噪方法的说明，本发明还提供一种降噪系统，其包括：发送端以及接收端，其中，发送端连续发出发射信号传送至大地系统，接收端接收经过大地系统传输后的响应信号。
其中，发送端连续发出的发射信号按照前述实施例1-实施例3等方法部分的陈述的要求。
在进一步可选的方案中，所述降噪系统还包括处理模块，其可以是设置于所述接收端内，或与所述接收端外并与所述接收端通信连接。譬如，接收端内部具有数据收发以及处理功能的处理器等。其利用接收的响应按照前述方法的方式进行降噪处理辨识出大地脉冲响应。或者是外置的电子终端设备，利用接收端反馈的数据进行辨识。
基于上述降噪方法的说明，本发明还提供一种可读存储介质，其存储了计算机程序，所述计算机程序被处理器调用以执行所述一种应用于可控源电磁法相关辨识中的降噪方法的步骤。
所述可读存储介质为计算机可读存储介质，其可以是前述任一实施例所述的控制器的内部存储单元，例如控制器的硬盘或内存。所述可读存储介质也可以是所述控制器的外部存储设备，例如所述控制器上配备的插接式硬盘，智能存储卡(Smart Media Card,SMC)，安全数字(Secure Digital,SD)卡，闪存卡(Flash Card)等。进一步地，所述可读存储介质还可以既包括所述控制器的内部存储单元也包括外部存储设备。所述可读存储介质用于存储所述计算机程序以及所述控制器所需的其他程序和数据。所述可读存储介质还可以用于暂时地存储已经输出或者将要输出的数据。
需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。
Claims (6)
1.一种应用于可控源电磁法相关辨识中的降噪方法，其特征在于：包括如下步骤：
设定发射信号的时长T与工频噪声周期Tc的关系；
连续发出发射信号传送至大地系统，并基于T与Tc的关系确定等效发射信号，其中，基于等效发射信号对应的等效接收信号中不包含工频噪声；
若发射信号的时长T为工频噪声周期Tc的整数倍，等效发射信号为连续发射的相邻两组发射信号相减得到；若发射信号的时长T为工频噪声周期Tc的整数倍加二分之一个Tc，等效发射信号为连续发射的相邻两组发射信号相加得到；
所述发射信号的时长T等于发射信号的有效激励时间与零激励时间之和；
若发射信号为m序列信号，发射信号的有效激励时间以及零激励时间对应的补零的个数如下所示：
Tv＝(Np*Lm*Nspb)/fs
式中，Npz为零激励时间对应的补零的个数，Tpz为零激励时间，Tv为发射信号的有效激励时间，Np为m序列的周期数，Lm为一个m序列周期的码元数，Nspb为每个码元的采样数，fs为采样频率，K为正整数；
其中，case1表示发射信号的时长T为工频噪声周期Tc的K倍，case2表示发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc；
还包括利用等效发射信号和等效接收信号辨识大地脉冲响应，其中，等效发射信号需满足：
case1：当发射信号的时长T为工频噪声周期Tc的K倍时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝-xm(t)，对应得到的等效发射信号为m序列2xm(t)；
case2：当发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝xm(t)，对应得到的等效发射信号为m序列2xm(t)，其中，xm(t)为m序列。
2.根据权利要求1所述的方法，其特征在于：所述辨识大地脉冲响应是利用m序列自相关特点采用维纳滤波反卷积法计算得到，其中，大地脉冲响应ht如下所示：
其中，为等效接收信号与等效发射信号的互相关函数，/>为等效发射信号的自相关函数；符号“*”表示卷积；采用反卷积求解大地脉冲响应ht存在：
其中，Aτ是等效发射信号st的归一化自相关值，Bτ是等效发射信号st和等效接收信号Vt的归一化互相关值，τ∈{0,1...,n}，n为采样点数，其分别表示如下：
其中，是等效接收信号Vt的自相关函数在零延时处的相关值，等效发射信号st和等效接收信号Vt分别表示如下：
y1(t)和y2(t)为连续发射的相邻两组发射信号x1(t)、x2(t)对应的响应信号。
3.根据权利要求1所述的降噪方法，其特征在于：case1：当发射信号的时长T为工频噪声周期Tc的K倍时，等效发射信号的功率为设计的相邻两段发射功率总和；
case2：当发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc时，等效发射信号的功率等于设计的发射功率。
4.根据权利要求1所述的降噪方法，其特征在于：若信号的发射端或响应信号的接收端存在固有误差，选择发射信号的时长T为工频噪声周期Tc的整数倍，等效发射信号为连续发射的相邻两组发射信号相减得到。
5.一种应用于可控源电磁法相关辨识中的降噪系统，其特征在于：包括发送端以及接收端，其中，发送端连续发出发射信号传送至大地系统，接收端接收经过大地系统传输后的响应信号；
若发射信号的时长T为工频噪声周期Tc的整数倍，连续发射的相邻两组发射信号相减得到等效发射信号；若发射信号的时长T为工频噪声周期Tc的整数倍加二分之一个Tc，连续发射的相邻两组发射信号相加得到等效发射信号，其中，基于等效发射信号对应的等效接收信号中不包含工频噪声；
所述发射信号的时长T等于发射信号的有效激励时间与零激励时间之和；
若发射信号为m序列信号，发射信号的有效激励时间以及零激励时间对应的补零的个数如下所示：
Tv＝(Np*Lm*Nspb)/fs
式中，Npz为零激励时间对应的补零的个数，Tpz为零激励时间，Tv为发射信号的有效激励时间，Np为m序列的周期数，Lm为一个m序列周期的码元数，Nspb为每个码元的采样数，fs为采样频率，K为正整数；
其中，case1表示发射信号的时长T为工频噪声周期Tc的K倍，case2表示发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc；
还包括处理模块，用于利用等效发射信号和等效接收信号辨识大地脉冲响应，其中，等效发射信号需满足：
case1：当发射信号的时长T为工频噪声周期Tc的K倍时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝-xm(t)，对应得到的等效发射信号为m序列2xm(t)；
case2：当发射信号的时长T为工频噪声周期Tc的K倍加二分之一个Tc时，连续发射的相邻两组发射信号x1(t)、x2(t)分别表示为：x1(t)＝xm(t)、x2(t)＝xm(t)，对应得到的等效发射信号为m序列2xm(t)，其中，xm(t)为m序列；
所述处理模块设置于所述接收端内，或与所述接收端外并与所述接收端通信连接。
6.一种可读存储介质，其特征在于：存储了计算机程序，所述计算机程序被处理器调用以执行权利要求1-4任一项所述方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110038837.5A CN112765546B (zh) | 2021-01-12 | 2021-01-12 | 一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110038837.5A CN112765546B (zh) | 2021-01-12 | 2021-01-12 | 一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112765546A CN112765546A (zh) | 2021-05-07 |
| CN112765546B true CN112765546B (zh) | 2024-05-14 |
Family
ID=75701694
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110038837.5A Active CN112765546B (zh) | 2021-01-12 | 2021-01-12 | 一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112765546B (zh) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106059684A (zh) * | 2016-07-14 | 2016-10-26 | 中国地质大学（北京） | 一种基于相关辨识的噪声处理系统及其方法 |
| CN107582068A (zh) * | 2017-10-12 | 2018-01-16 | 南阳师范学院 | 一种基于周期自相干原理的畸变产物耳声发射探测方案 |
| CN109001828A (zh) * | 2018-08-08 | 2018-12-14 | 西安石油大学 | 一种地下目标体的识别方法及装置、存储介质 |
| CN110398778A (zh) * | 2019-07-10 | 2019-11-01 | 吉林大学 | 一种基于等效采样的航空电磁浅层数据相关辨识方法 |
| CN111767887A (zh) * | 2020-07-08 | 2020-10-13 | 吉林大学 | 基于小波分解与ime频率估计的瞬变电磁数据处理方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120179372A1 (en) * | 2010-07-22 | 2012-07-12 | Alexander Edward Kalish | Collecting Control Source Electromagnetic Signals |
-
2021
- 2021-01-12 CN CN202110038837.5A patent/CN112765546B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106059684A (zh) * | 2016-07-14 | 2016-10-26 | 中国地质大学（北京） | 一种基于相关辨识的噪声处理系统及其方法 |
| CN107582068A (zh) * | 2017-10-12 | 2018-01-16 | 南阳师范学院 | 一种基于周期自相干原理的畸变产物耳声发射探测方案 |
| CN109001828A (zh) * | 2018-08-08 | 2018-12-14 | 西安石油大学 | 一种地下目标体的识别方法及装置、存储介质 |
| CN110398778A (zh) * | 2019-07-10 | 2019-11-01 | 吉林大学 | 一种基于等效采样的航空电磁浅层数据相关辨识方法 |
| CN111767887A (zh) * | 2020-07-08 | 2020-10-13 | 吉林大学 | 基于小波分解与ime频率估计的瞬变电磁数据处理方法 |
Non-Patent Citations (5)
| Title |
|---|
| "Controlled-source electromagnetic data processing based on gray system theory and robust estimation ";Dan Mo等;《Applied Geophysics》;20180214;第14卷(第4期);第570-580页 * |
| "电磁法相关辨识抗干扰机理研究与数值模拟";智庆全 等;《地球物理学进展》;20180123;第33卷(第4期);第1755-1766页 * |
| "相关辨识的算法优化";甄晓丹;《中国优秀硕士学位论文全文数据库基础科学辑》;20180815;第A003-2页 * |
| 基于相关估计辨识工频噪声中的微弱信号;曹才开 等;《电气电子教学学报》;20030325;第25卷(第2期);第59-61+64页 * |
| 杨智 ; 等."基于双m序列的可控源电磁法强周期噪声对消方法研究".《中国地球科学联合学术年会2021》.2021,第1851页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112765546A (zh) | 2021-05-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yadav et al. | Electrocardiogram signal denoising using non‐local wavelet transform domain filtering | |
| Wu et al. | Elimination of cross‐terms in the Wigner–Ville distribution of multi‐component LFM signals | |
| CN110690931B (zh) | 一种基于多小波基联合的数字信号自适应码率估计方法和装置 | |
| CN113972929B (zh) | 一种高动态多普勒下扩频信号的捕获方法 | |
| Yoo et al. | Image restoration by blind‐Wiener filter | |
| Orović et al. | Improved higher order robust distributions based on compressive sensing reconstruction | |
| CN112765546B (zh) | 一种应用于可控源电磁法相关辨识中的降噪方法、系统及可读存储介质 | |
| CN117692074B (zh) | 一种适用于非稳态水声目标信号的低频混叠噪声抑制方法 | |
| CN108810923B (zh) | Afh干扰频点的判断方法及装置、计算机可读存储介质、接收机 | |
| CN109765608B (zh) | 一种基于联合字典的煤层巷道锚杆振动噪声压制方法 | |
| Larsen et al. | Suppression of very low frequency radio noise in transient electromagnetic data with semi-tapered gates | |
| CN107704831A (zh) | 一种基于奇异值分解中值法的瓦斯浓度数据降噪方法 | |
| CN112630824B (zh) | 一种地震成像中的离散点扩散函数生成方法及系统 | |
| CN110632563B (zh) | 一种基于短时傅里叶变换的脉内频率编码信号参数测量方法 | |
| CN116679325A (zh) | 一种新的gnss快扫频干扰抑制方法 | |
| CN116953631A (zh) | 全极化雷达发射波形优化方法、计算机装置和存储介质 | |
| CN113341463B (zh) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 | |
| CN110417696A (zh) | 一种基于随机共振的2ask信号解调方法 | |
| CN114611542A (zh) | 信号降噪处理方法及通信装置 | |
| Zhao et al. | Robust transcale state estimation for multiresolution discrete‐time systems based on wavelet transform | |
| Tadokoro et al. | A simple optimum nonlinear filter for stochastic-resonance-based signal detection | |
| Treichler | γ-LMS and its use in a noise-compensating adaptive spectral analysis technique | |
| CN111325959A (zh) | 一种无线无源型传感器信号检测电路和方法 | |
| CN113504526B (zh) | 基于mimo雷达的目标检测方法、装置、电子设备及存储介质 | |
| CN111384983B (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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |