CN112285793B - A kind of magnetotelluric denoising method and system - Google Patents
A kind of magnetotelluric denoising method and system Download PDFInfo
- Publication number
- CN112285793B CN112285793B CN202011309115.0A CN202011309115A CN112285793B CN 112285793 B CN112285793 B CN 112285793B CN 202011309115 A CN202011309115 A CN 202011309115A CN 112285793 B CN112285793 B CN 112285793B
- Authority
- CN
- China
- Prior art keywords
- data segment
- matrix
- denoised
- sounding
- magnetotelluric
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 51
- 239000011159 matrix material Substances 0.000 claims abstract description 260
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 29
- 230000010349 pulsation Effects 0.000 claims abstract 19
- 239000013598 vector Substances 0.000 claims description 99
- 238000012545 processing Methods 0.000 claims description 28
- 238000012216 screening Methods 0.000 claims description 17
- 230000011218 segmentation Effects 0.000 claims description 15
- 238000005259 measurement Methods 0.000 claims description 14
- 238000010276 construction Methods 0.000 claims description 8
- 238000005516 engineering process Methods 0.000 description 4
- 238000000926 separation method Methods 0.000 description 4
- 239000006185 dispersion Substances 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000005672 electromagnetic field Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000005065 mining Methods 0.000 description 2
- 230000010287 polarization Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 238000007430 reference method Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 230000001629 suppression Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明涉及一种大地电磁去噪方法及系统，对大地电磁脉动信号进行重叠分段，确定每个测深数据段的近似熵，将近似熵大于近似熵阈值的测深数据段确定为有效测深数据段，将其他测深数据段确定为待去噪测深数据段，选用近似熵对数据段进行筛选，防止了对有效测深数据段的“过处理”；对待去噪测深数据段进行奇异值分解，将分解后的奇异值矩阵、左和右奇异矩阵相乘，获得多个分量矩阵，构造相应的分量数据段，确定每个分量数据段的近似熵，将筛选的近似熵大于测深数据段的近似熵的分量数据段求和，获得去噪后的测深数据段，实现了去噪的目的；最后将去噪后的测深数据段和有效测深数据段重构，获得重构去噪大地电磁脉动信号，提高了数据去噪的准确性。
The invention relates to a magnetotelluric denoising method and system. The magnetotelluric pulsation signal is overlapped and segmented, the approximate entropy of each sounding data segment is determined, and the sounding data segment whose approximate entropy is greater than the approximate entropy threshold is determined as an effective sounding data segment. For the depth data segment, other sounding data segments are determined as the sounding data segments to be denoised, and approximate entropy is used to filter the data segments, which prevents the "overprocessing" of the effective sounding data segments; Perform singular value decomposition, multiply the decomposed singular value matrix, left and right singular matrix, obtain multiple component matrices, construct corresponding component data segments, determine the approximate entropy of each component data segment, and filter the approximate entropy greater than The approximate entropy component data segments of the sounding data segment are summed to obtain the denoised sounding data segment, which realizes the purpose of denoising. It can denoise the magnetotelluric pulsation signal and improve the accuracy of data denoising.
Description
Technical Field
The invention relates to the field of data processing, in particular to a magnetotelluric denoising method and a magnetotelluric denoising system.
Background
A Magnetotelluric sounding Method (MT) observes a component disturbance value generated by a natural electromagnetic field signal through an underground medium on the ground surface so as to acquire an electrical structure of the underground medium. The field source is a natural electromagnetic field, effective signals originate from interaction of lightning discharge and solar wind with an earth magnetic layer, and compared with an electromagnetic sounding method of an artificial field source, the MT has the advantages of rich frequency band information, large detection depth and the like, so that the MT is widely applied to a plurality of fields of oil and gas field general survey exploration, earthquake prediction, mine exploration and the like. However, the natural electromagnetic field signal is weak, and is easily interfered by various electromagnetic noises during the measurement, especially in the strong interference environment such as a mining area. The effective signal has the obvious inconsistent property with noise interference, the effective signal has the characteristics of randomness, non-stationarity, non-Gaussian and the like, the interference signal usually has a certain polarization direction and strong signal energy, and when the strong interference source is close to the measuring point, the interference signal can completely cover the effective signal. Thus, the data quality is very poor, and serious deviation occurs in estimating the impedance, thereby causing a near source interference phenomenon. Therefore, achieving effective suppression of noise is a prerequisite for the MT method to be able to obtain an accurate electrical structure.
The existing methods for removing magnetotelluric noise can be mainly classified into the following three categories: the time domain processing method comprises the following steps: the method mainly identifies and removes noise forms obviously existing in a time sequence, and typically represents a method such as a form filtering method, a subspace enhancement algorithm, matching pursuit, compressed sensing and the like. For example, a morphological filtering method using rectangular structural elements only has the best processing effect on square waves, but cannot effectively remove large-scale noise and abrupt noise in other forms such as charge-discharge triangular waves. Therefore, at present, a plurality of time domain methods are often combined to perform denoising processing, but the corresponding processing time is increased; the time-frequency conversion method comprises the following steps: the time series is converted into other domains (such as frequency domain, wavelet domain, etc.), noise is removed by screening the spectral information of the signal, which typically represents methods such as wavelet analysis, Empirical mode decomposition (Empirical mode decomposition), and Variational mode decomposition (Variational mode decomposition), but has the problems of mother wavelet function and decomposition layer number selection, modal aliasing, preset parameter selection, etc. Thirdly, the frequency domain processing method comprises the following steps: the method comprises a least square method, Robust estimation, a far reference method and the like, and is mainly associated with impedance estimation, for example, the far reference method replaces self-power spectrum data of a measuring point by cross-power spectrum data between the reference point and the measuring point to perform impedance calculation, and can remove noise irrelevant to the reference point in the measuring point data. Therefore, when the signal-to-noise ratio of the acquired data is high or the form of the existing noise is regular, the existing magnetotelluric denoising technology can achieve a good processing effect, but in actual observation, the noise form is often complex and changeable, and particularly in a mine collection area, the requirement of high signal-to-noise ratio is difficult to meet.
Based on the difference between the properties of the magnetotelluric effective signal and the interfering signal, it is proposed to remove the noise by using SVD (Singular Value Decomposition). Before SVD, a one-dimensional time sequence vector needs to be converted into a multi-dimensional matrix, the SVD extracts different singular values (arranged from large to small and can be understood as field source values to a certain extent) through the overall characteristics of the multi-dimensional matrix and the local characteristics of each row of vectors, calculates each decomposition component of the time sequence vector through the singular values, each component is an electromagnetic wave signal excited by different field sources, and achieves the purpose of denoising by judging whether the component is an effective signal or not. Therefore, the method is mainly related to the number of decomposed components when removing the noise (when the number of the components is large, the more detailed signal decomposition is, the effective signal and the noise can be completely separated, but the processing time is correspondingly increased; when the number of the components is small, the processing speed is high, but the effective signal can be lost), has no requirement on the signal-to-noise ratio, and can be applied to data denoising processing of a mining area. And the method does not carry out denoising processing according to the noise form, but identifies and removes the interference (singular value) field source. The electromagnetic wave form excited by the interference source is complex and changeable, so that the method can simultaneously process various types of noise and is not limited by the noise form. However, when the data segment is a pure signal, the denoising processing of the screening component is still performed, and the problem of "over processing" of the effective signal occurs.
Disclosure of Invention
The invention aims to provide a magnetotelluric denoising method and a magnetotelluric denoising system, which are used for reducing the over-processing of data and improving the data denoising accuracy.
In order to achieve the purpose, the invention provides the following scheme:
a magnetotelluric denoising method, the method comprising:
measuring magnetotelluric pulse signals of different depths of the earth by using a magnetotelluric sounding method;
carrying out overlapping segmentation on the magnetotelluric pulse signals to obtain a plurality of sounding data segments;
determining the approximate entropy of each sounding data segment based on the approximate entropy theory;
determining the depth measuring data segment with the approximate entropy smaller than or equal to the approximate entropy threshold value as a depth measuring data segment to be denoised, and determining the depth measuring data segment with the approximate entropy larger than the approximate entropy threshold value as an effective depth measuring data segment;
constructing a matrix to be denoised according to the sounding data segment to be denoised;
performing singular value decomposition on the matrix to be denoised to obtain a singular value matrix, a left singular matrix and a right singular matrix of the matrix to be denoised;
multiplying each diagonal element of the singular value matrix, the column vector corresponding to the left singular matrix and the row vector corresponding to the right singular matrix to obtain a plurality of field source component matrixes of the sounding data segment to be denoised;
converting each field source component matrix into a component data segment;
determining the approximate entropy of each component data segment based on the approximate entropy theory;
screening the component data sections with the approximate entropy larger than that of the depth measuring data section to be denoised, and summing the screened component data sections to obtain the denoised depth measuring data section;
and reconstructing the denoised sounding data segment and the effective sounding data segment to obtain a reconstructed denoised magnetotelluric pulse signal.
Optionally, the performing overlapping segmentation on the magnetotelluric pulse signal to obtain a plurality of sounding data segments specifically includes:
for the magnetotelluric pulse signal x1,x2,…xwPerforming overlapping segmentation processing to obtain each depth measurement data segment Xv＝{x(v-1)*N+1,x(v-1)*N+2,…,x(v-1)*N+N}；
Wherein x is1,x2,…xwRespectively the 1 st, 2 nd and w th magnetotelluric pulse signals, w is the length of the magnetotelluric pulse signal, x(v-1)*N+1、x(v-1)*N+2、x(v-1)*N+NThe data are respectively the 1 st, 2 nd and nth data of the v-th sounding data segment, N is the sounding data segment length, v is 1,2, …, l is the number of sounding data segments, l is (w- α N)/N (1- α), and α is the data overlap ratio.
Optionally, the determining the approximate entropy of each depth measurement data segment based on the approximate entropy theory specifically includes:
initializing the embedding dimension of the Henkel matrix as h;
constructing a Hankel matrix of the embedding dimension for each sounding data segment; each column in the hankel matrix is a mode vector;
using formulas
Determining the maximum distance between the p mode vector in the Henkel matrix and the q mode vector of the Henkel matrix;
counting the number of the p-th mode vector in the Hankel matrix, wherein the maximum distance between the p-th mode vector and all the mode vectors of the Hankel matrix is smaller than or equal to a mode fault tolerance threshold;
based on the approximate probability theory, according to the number, using formula Cp h(r)＝Np h(r)/(N-h +1), determining an approximate probability value of the p-th mode vector in the hankel matrix;
using a formula based on the approximate probability value
Determining an approximate probability average value of each sounding data segment under the h embedding dimension;
increasing the embedding dimension h by 1, returning to the step of constructing a Hankel matrix with the embedding dimension h for each depth measuring data segment, and determining the approximate probability average value of each depth measuring data segment under the embedding dimension h + 1;
according to the approximate probability average value of each sounding data segment in the h embedding dimension and the approximate probability average value of each sounding data segment in the h +1 embedding dimension, a formula is utilized
Determining an approximate entropy of each sounding data segment;
wherein d ismax[y(p),y(q)]Is the maximum distance between the p-th mode vector in the Hankel matrix and the q-th mode vector in the Hankel matrix, y (p) is the p-th mode vector in the Hankel matrix, y (q) is the q-th mode vector in the Hankel matrix, and x (p + delta) is the Hankel matrixThe delta-th element of the p-th mode vector in the matrix, x (q + delta) is the delta-th element of the q-th mode vector in the hank matrix, delta is 1,2, …, h, Cp h(r) is the approximate probability value of the p-th pattern vector in the Hankel matrix, Np h(r) is the number of maximum distances less than or equal to the mode fault tolerance threshold, r is the mode fault tolerance threshold, N is the sounding data segment length,
is an approximate probability average of the sounding data segment in the h embedding dimension,
is the approximate probability average, Ae, of the sounding data segment in the h +1 embedding dimensioniIs the approximate entropy of the ith sounding data segment.
Optionally, constructing a matrix to be denoised according to the sounding data segment to be denoised specifically includes:
according to the preset matrix dimension, the a-th sounding data segment X to be denoiseda＝{x(a-1)*N+1,x(a-1)*N+2,…,x(a-1)*N+NEqually dividing into a plurality of sections;
constructing a matrix to be denoised according to each section after the equal division
Wherein, XaFor the a-th sounding data segment to be denoised, x(a-1)*N+1、x(a-1)*N+2、x(a-1)*N+n、x(a-1)*N+1+n、x(a-1)*N+2n、xa*N+1-nAnd x(a-1)*N+NRespectively representing the 1 st, 2 nd, nth, N +1 st, 2 nth, N-N +1 th and nth data of the depth measuring data segment to be denoised, wherein N is the length of the depth measuring data segment, N is the length of a preset matrix, N is N/m, and m is the dimension of the preset matrix.
Optionally, the performing singular value decomposition on the matrix to be denoised to obtain a singular value matrix, a left singular matrix, and a right singular matrix of the matrix to be denoised, and then further includes:
determining the ratio of each singular value in the singular value matrix to the sum of all singular values;
drawing a singular value proportion broken line graph by taking the number of singular values in the singular value matrix as an abscissa and the ratio as an ordinate;
judging whether the vertical coordinates corresponding to the tail branches of the singular value proportion broken line graph are all smaller than a ratio threshold value or not, and obtaining a judgment result;
if the judgment result shows that the matrix to be denoised is the zero matrix, outputting a singular value matrix, a left singular matrix and a right singular matrix of the matrix to be denoised;
if the judgment result shows no, updating the preset matrix dimension, and returning to the step of' according to the preset matrix dimension, enabling the a-th depth sounding data segment X to be denoiseda＝{x(a-1)*N+1,x(a-1)*N+2,…,x(a-1)*N+NAnd equally dividing into multiple sections.
Optionally, reconstructing the denoised sounding data segment and the effective sounding data segment to obtain a reconstructed denoised magnetotelluric pulse signal, specifically including:
reconstructing the denoised sounding data segment and the effective sounding data segment by using a formula w ═ lN- (l-1) · α N to obtain a reconstructed denoised magnetotelluric pulse signal;
wherein w is the length of the magnetotelluric pulse signal, l is the number of the sounding data segments, alpha is the data overlapping rate, and N is the length of the sounding data segments.
A magnetotelluric denoising system, the system comprising:
the magnetotelluric pulse signal measuring module is used for measuring magnetotelluric pulse signals of different depths of the earth by utilizing a magnetotelluric sounding method;
a sounding data segment obtaining module, configured to perform overlapping segmentation on the magnetotelluric pulse signals to obtain multiple sounding data segments;
the approximate entropy determining module of the sounding data segments is used for determining the approximate entropy of each sounding data segment based on the approximate entropy theory;
the depth measuring data segment determining module is used for determining a depth measuring data segment with approximate entropy smaller than or equal to the approximate entropy threshold value as a depth measuring data segment to be denoised and determining a depth measuring data segment with approximate entropy larger than the approximate entropy threshold value as an effective depth measuring data segment;
the matrix construction module to be denoised is used for constructing a matrix to be denoised according to the sounding data segment to be denoised;
the singular value decomposition module is used for performing singular value decomposition on the matrix to be denoised to obtain a singular value matrix, a left singular matrix and a right singular matrix of the matrix to be denoised;
a field source component matrix obtaining module, configured to multiply each diagonal element of the singular value matrix, a column vector corresponding to the left singular matrix, and a row vector corresponding to the right singular matrix, to obtain a plurality of field source component matrices of the sounding data segment to be denoised;
a component data segment obtaining module for converting each field source component matrix into a component data segment;
the approximate entropy determining module of the component data segments is used for determining the approximate entropy of each component data segment based on the approximate entropy theory;
the denoised depth measurement data section obtaining module is used for screening the component data sections of which the approximate entropy is larger than that of the depth measurement data section to be denoised, and summing the screened component data sections to obtain the denoised depth measurement data section;
and the reconstruction denoising magnetotelluric pulse signal obtaining module is used for reconstructing the denoised sounding data segment and the effective sounding data segment to obtain a reconstruction denoising magnetotelluric pulse signal.
Optionally, the depth measurement data segment obtaining module specifically includes:
a sounding data segment obtaining submodule for obtaining the magnetotelluric pulse signal x1,x2,…xwPerforming overlapping segmentation processing to obtain each depth measurement data segment Xv＝{x(v-1)*N+1,x(v-1)*N+2,…,x(v-1)*N+N}；
Wherein x is1,x2,…xwRespectively the 1 st, 2 nd and w th magnetotelluric pulse signals, w is the length of the magnetotelluric pulse signal, x(v-1)*N+1、x(v-1)*N+2、x(v-1)*N+NThe data are respectively the 1 st, 2 nd and nth data of the v-th sounding data segment, N is the sounding data segment length, v is 1,2, …, l is the number of sounding data segments, l is (w- α N)/N (1- α), and α is the data overlap ratio.
Optionally, the approximate entropy determining module of the sounding data segment specifically includes:
the embedded dimension initialization submodule is used for initializing the embedded dimension of the Hankel matrix to be h;
the Hankel matrix construction submodule is used for respectively constructing a Hankel matrix with an embedded dimension h for each sounding data segment; each column in the hankel matrix is a mode vector;
maximum distance determination submodule for utilizing a formula
Determining the maximum distance between the p mode vector in the Henkel matrix and the q mode vector of the Henkel matrix;
the number counting submodule is used for counting the number of the p-th mode vector in the Hankel matrix, wherein the maximum distance between the p-th mode vector and all the mode vectors of the Hankel matrix is smaller than or equal to a mode fault tolerance threshold;
an approximate probability value determining submodule for utilizing a formula C according to the number based on the approximate probability theoryp h(r)＝Np h(r)/(N-h +1), determining an approximate probability value of the p-th mode vector in the hankel matrix;
an approximate probability average determination submodule for utilizing a formula according to the approximate probability value
Determining an approximate probability average value of each sounding data segment under the h embedding dimension;
the circulation submodule is used for increasing the embedding dimension h by 1, returning to the step of 'respectively constructing a Hankel matrix with the embedding dimension h for each sounding data segment', and determining the approximate probability average value of each sounding data segment under the embedding dimension h + 1;
an approximate entropy determination submodule of the sounding data segment for utilizing a formula according to the approximate probability average value of each sounding data segment in the h embedding dimension and the approximate probability average value of each sounding data segment in the h +1 embedding dimension
Determining an approximate entropy of each sounding data segment;
wherein d ismax[y(p),y(q)]Is the maximum distance between the p-th mode vector in the Heckel matrix and the q-th mode vector in the Heckel matrix, y (p) is the p-th mode vector in the Heckel matrix, y (q) is the q-th mode vector in the Heckel matrix, x (p + delta) is the delta-th element of the p-th mode vector in the Heckel matrix, x (q + delta) is the delta-th element of the q-th mode vector in the Heckel matrix, delta is 1,2, …, h, Cp h(r) is the approximate probability value of the p-th pattern vector in the Hankel matrix, Np h(r) is the number of maximum distances less than or equal to the mode fault tolerance threshold, r is the mode fault tolerance threshold, N is the sounding data segment length,
is an approximate probability average of the sounding data segment in the h embedding dimension,
is the approximate probability average, Ae, of the sounding data segment in the h +1 embedding dimensioniIs the approximate entropy of the ith sounding data segment.
Optionally, the matrix construction module to be denoised specifically includes:
an equal division module used for dividing the a-th depth sounding data segment X to be denoised according to the preset matrix dimensiona＝{x(a-1)*N+1,x(a-1)*N+2,…,x(a-1)*N+NEqually dividing into a plurality of sections;
a sub-module for constructing the matrix to be denoised according to each segment after being equally divided
Wherein, XaFor the a-th sounding data segment to be denoised, x(a-1)*N+1、x(a-1)*N+2、x(a-1)*N+n、x(a-1)*N+1+n、x(a-1)*N+2n、xa*N+1-nAnd x(a-1)*N+NRespectively representing the 1 st, 2 nd, nth, N +1 st, 2 nth, N-N +1 th and nth data of the depth measuring data segment to be denoised, wherein N is the length of the depth measuring data segment, N is the length of a preset matrix, N is N/m, and m is the dimension of the preset matrix.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
the invention discloses a magnetotelluric denoising method and a magnetotelluric denoising system, wherein magnetotelluric pulse signals are subjected to overlapping segmentation to obtain a plurality of sounding data sections, the approximate entropy of each sounding data section is determined, the sounding data section with the approximate entropy smaller than or equal to the approximate entropy threshold is determined as a sounding data section to be denoised, the sounding data section with the approximate entropy larger than the approximate entropy threshold is determined as an effective sounding data section, and the data sections are screened by selecting the approximate entropy before performing signal-noise separation by adopting SVD (singular value decomposition), so that the effective sounding data section is prevented from being over-processed; carrying out singular value decomposition on the sounding data segment to be denoised, multiplying a singular value matrix, a left singular matrix and a right singular matrix after the singular value decomposition to obtain a plurality of component matrices of the sounding data segment to be denoised, converting each component matrix into a component data segment, determining the approximate entropy of each component data segment, screening the component data segments of which the approximate entropy is greater than that of the sounding data segment to be denoised, summing the screened component data segments to obtain the denoised sounding data segment, and screening the decomposed components by adopting the approximate entropy to realize the purpose of denoising; and finally reconstructing the denoised depth sounding data segment and the effective depth sounding data segment to obtain a reconstructed denoised magnetotelluric pulse signal, so that the accuracy of data denoising is improved.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed to be used in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without inventive exercise.
FIG. 1 is a flow chart of a magnetotelluric denoising method provided by the present invention;
fig. 2 is a schematic diagram of a magnetotelluric denoising method provided by the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
The invention aims to provide a magnetotelluric denoising method and a magnetotelluric denoising system, which are used for reducing the over-processing of data and improving the data denoising accuracy.
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.
The invention provides a magnetotelluric denoising method in order to reduce the loss of effective signals and combine approximate entropy screening and SVD decomposition denoising, as shown in fig. 1 and fig. 2, the method comprises the following steps:
and S101, measuring magnetotelluric pulse signals of different depths of the earth by using a magnetotelluric sounding method.
And S102, carrying out overlapping segmentation on the magnetotelluric pulse signals to obtain a plurality of sounding data segments.
S103, determining the approximate entropy of each sounding data segment based on the approximate entropy theory.
S104, determining the depth measuring data segment with the approximate entropy smaller than or equal to the approximate entropy threshold value as a depth measuring data segment to be denoised, and determining the depth measuring data segment with the approximate entropy larger than the approximate entropy threshold value as an effective depth measuring data segment. Preferably, the approximate entropy threshold is Ae0，Ae0Typically between 0.5 and 1.
And S105, constructing a matrix to be denoised according to the sounding data segment to be denoised.
And S106, performing singular value decomposition on the matrix to be denoised to obtain a singular value matrix, a left singular matrix and a right singular matrix of the matrix to be denoised.
And S107, multiplying each diagonal element of the singular value matrix, the column vector corresponding to the left singular matrix and the row vector corresponding to the right singular matrix to obtain a plurality of field source component matrixes of the sounding data segment to be denoised.
And S108, converting each field source component matrix into a component data segment.
S109, based on the approximate entropy theory, the approximate entropy of each component data segment is determined.
S110, screening the component data sections of which the approximate entropy is larger than that of the depth measuring data section to be denoised, and summing the screened component data sections to obtain the denoised depth measuring data section.
And S111, reconstructing the denoised sounding data segment and the effective sounding data segment to obtain a reconstructed denoised magnetotelluric pulse signal.
The specific process is as follows:
step S102, carrying out overlapping segmentation on the magnetotelluric pulse signals to obtain a plurality of sounding data segments, which specifically comprises the following steps:
to magnetotelluric pulsating signal x1,x2,…xwPerforming overlapping segmentation processing to obtain each depth measurement data segment Xv＝{x(v-1)*N+1,x(v-1)*N+2,…,x(v-1)*N+N}。
Wherein x is1,x2,…xwRespectively the 1 st, 2 nd and w th magnetotelluric pulsesSignal, w is the length of the magnetotelluric pulse signal, x(v-1)*N+1、x(v-1)*N+2、x(v-1)*N+NThe data are respectively the 1 st, 2 nd and nth data of the v-th sounding data segment, N is the sounding data segment length, v is 1,2, …, l is the number of sounding data segments, l is (w- α N)/N (1- α), and α is the data overlap ratio.
Because the multidimensional matrix of the processing object decomposed by SVD has larger requirement on the running memory when the processing object is processed by computer software, the step firstly divides the one-dimensional vector into a plurality of data segments, and each data segment is constructed into the multidimensional matrix, thereby shortening the length of the processing object, and ensuring the feasibility of running, the continuity of data and the accuracy of screening the data segments.
S103, determining the approximate entropy of each sounding data segment based on the approximate entropy theory, and specifically comprising the following steps:
initializing an embedding dimension h of a Hankel matrix (Hankel matrix);
constructing a Hankel matrix with embedded dimensions for each sounding data segment; each column in the hankel matrix is a mode vector;
using formulas
Determining the maximum distance between the p-th mode vector in the Hankel matrix and the q-th mode vector of the Hankel matrix;
counting the number of the p-th mode vector in the Hankel matrix, wherein the maximum distance between the p-th mode vector and all the mode vectors of the Hankel matrix is smaller than or equal to the mode fault-tolerant threshold;
based on approximate probability theory, according to the number, using formula Cp h(r)＝Np h(r)/(N-h +1), determining the approximate probability value of the p-th mode vector in the hankel matrix;
using a formula based on the approximate probability value
Determining an approximate probability average of each sounding data segment in the h-embedding dimension；
Increasing the embedding dimension h by 1, returning to the step of constructing a Hankel matrix with the embedding dimension h for each depth measuring data segment, and determining the approximate probability average value of each depth measuring data segment under the embedding dimension h + 1;
according to the approximate probability average value of each sounding data segment in the h embedding dimension and the approximate probability average value of each sounding data segment in the h +1 embedding dimension, a formula is utilized
Determining an approximate entropy of each sounding data segment;
wherein d ismax[y(p),y(q)]Is the maximum distance between the p-th mode vector in the Heckel matrix and the q-th mode vector in the Heckel matrix, y (p) is the p-th mode vector in the Heckel matrix, y (q) is the q-th mode vector in the Heckel matrix, x (p + delta) is the delta-th element of the p-th mode vector in the Heckel matrix, x (q + delta) is the delta-th element of the q-th mode vector in the Heckel matrix, delta is 1,2, …, h, Cp h(r) is the approximate probability value of the p-th pattern vector in the Hankel matrix, Np h(r) is the number of maximum distances less than or equal to the mode fault tolerance threshold, r is the mode fault tolerance threshold, N is the sounding data segment length,
is an approximate probability average of the sounding data segment in the h embedding dimension,
is the approximate probability average, Ae, of the sounding data segment in the h +1 embedding dimensioniIs the approximate entropy of the ith sounding data segment.
The approximate entropy is adopted to carry out pre-screening on each segment of data before signal-noise separation, and noise has stronger directivity and more concentrated dispersion degree compared with effective signals, so that the data with better processing quality can be eliminated by the difference through the pre-screening by means of the approximate entropy, and the phenomenon that the effective signals are processed by denoising of the screening component when the data segment is the effective pure signals and the 'over-processing' of the screening component is still carried out is prevented.
Step S105, constructing a matrix to be denoised according to the sounding data segment to be denoised, which specifically comprises the following steps:
according to the preset matrix dimension, the a-th sounding data segment X to be denoiseda＝{x(a-1)*N+1,x(a-1)*N+2,…,x(a-1)*N+NEqually dividing into a plurality of sections;
constructing a matrix to be denoised according to each section after the equal division
Wherein, XaFor the a-th sounding data segment to be denoised, x(a-1)*N+1、x(a-1)*N+2、x(a-1)*N+n、x(a-1)*N+1+n、x(a-1)*N+2n、xa*N+1-nAnd x(a-1)*N+NRespectively representing the 1 st, 2 nd, nth, N +1 st, 2 nth, N-N +1 th and nth data of the depth measuring data segment to be denoised, wherein N is the length of the depth measuring data segment, N is the length of a preset matrix, N is N/m, and m is the dimension of the preset matrix.
Step S106, SVD decomposition: xa＝USVTWherein U ═ U1…Um]Is a left singular matrix of m x m orthogonal,
is a matrix of singular values, the principal diagonal is the singular value of the matrix,
is Xa TXaN × n right singular matrices. U shape1、UmRespectively 1 st and m-th elements, S, in the left singular matrix1、SmRespectively representing the 1 st and m-th elements, V, of the matrix of singular values1 T、V1 TRespectively representing the 1 st and nth elements in the right singular matrix.
The dimension and length of the matrix structure are selected very key during SVD decomposition, and the suppression effect on noise is better when the length of the matrix is less than or equal to the length of the noise through repeated tests.
Step S106, then further comprising:
determining the ratio of each singular value in the singular value matrix to the sum of all singular values; delta. in FIG. 2zFor the z-th singular value, delta, in the matrix of singular valuessumThe sum of singular values in a singular value matrix is obtained;
drawing a singular value proportion broken line graph by taking the number of singular values in the singular value matrix as an abscissa and the ratio as an ordinate;
judging whether the vertical coordinates corresponding to the tail branches of the singular value proportion broken line graph are all smaller than a ratio threshold value or not, and obtaining a judgment result; preferably, the threshold ratio is 0.001. The tail branch refers to the last part of the fold line in the singular value gravity fold line graph.
If the judgment result shows that the data segment is the deep sounding data segment to be denoised, outputting a singular value matrix, a left singular matrix and a right singular matrix of the deep sounding data segment to be denoised;
if the judgment result shows no, updating the preset matrix dimension, and returning to the step of' according to the preset matrix dimension, enabling the a-th sounding data segment X to be denoiseda＝{x(a-1)*N+1,x(a-1)*N+2,…,x(a-1)*N+NAnd equally dividing into multiple sections.
And judging whether the selection of the matrix dimension is correct or not by observing whether the tail branch of the line graph is less than 0.001 or not according to the singular value proportion line graph, and if the tail branch value of the line graph is more than 0.001, re-assigning the matrix dimension until the proportion value is less than 0.001.
Step S107, calculating all component matrixes according to the singular values and the left and right eigenvectors, connecting the component matrixes end to obtain signals transmitted by different field source information:
wherein, Xa 1、Xa mRespectively being the 1 st and the m-th component matrixes of the a-th sounding data segment to be denoised.
Step S108, the step of converting each field source component matrix in the field source matrix into a component data segment is the opposite step of "constructing a matrix to be denoised according to the sounding data segment to be denoised" in step S105.
In step S109, the method of determining the approximate entropy of each component data segment based on the approximate entropy theory is the same as that of step S103.
S110, the component data segment with the approximate entropy larger than that of the sounding data segment to be denoised is regarded as effective data to be reserved, and the component data segment with the approximate entropy smaller than that is regarded as noise data. And adding and summing the reserved effective components to realize the signal-noise separation of the sounding data section to be denoised. X in FIG. 2bRepresenting a denoised sounding data segment, Xb kThe k-th component matrix is shown, where k is 1,2, …, m.
Step S111, reconstructing the denoised sounding data segment and the effective sounding data segment to obtain a reconstructed denoised magnetotelluric pulse signal, which specifically comprises the following steps:
reconstructing the denoised depth sounding data segment and the effective depth sounding data segment by using a formula w ═ lN- (l-1) · α N to obtain a reconstructed denoised magnetotelluric pulse signal;
wherein w is the length of the magnetotelluric pulse signal, l is the number of the sounding data segments, alpha is the data overlapping rate, and N is the length of the sounding data segments.
The basic principle of the invention is as follows: selecting a data segment to be processed by using approximate entropy, decomposing a time sequence signal into a singular value matrix and left and right characteristic matrix vectors by using SVD (singular value decomposition) based on the property difference of a natural field source and an interference source, and then forming a plurality of component signals from different field source information by SVD inverse transformation, and calculating the approximate entropy of the component signals, wherein the approximate entropy of effective signals is larger (representing that the signal dispersion degree is larger); and the approximate entropy of the noise is smaller (with a certain polarization direction), so as to screen the components.
Therefore, the approximate entropy and the SVD are combined, so that the data over-processing can be reduced, and the component where the effective signal is located can be accurately selected to reconstruct the de-noised data segment. Under the conditions of low data quality and complex form noise of actually measured data, the magnetotelluric denoising technology based on approximate entropy and SVD decomposition can recover the change trend of an apparent resistivity-phase curve to obtain a continuous and smooth curve form.
According to the invention, the approximate entropy is adopted to pre-screen each segment of data before signal-noise separation, and noise has stronger directivity and more concentrated dispersion degree compared with effective signals, so that the data with better quality can be prevented from being further processed by pre-screening through the difference; and a plurality of singular values are extracted through the difference between the signal amplitudes, each singular value can calculate signal components containing different information, and the purpose of denoising is achieved through screening the components, so that the technology is irrelevant to the signal-to-noise ratio of data and is relevant to the number of the extracted singular values (the more the singular values are, the more the components are, the slower the processing speed is; the less the singular values are, the less the components are, effective signals can be lost); similarly, the invention does not carry out denoising according to the form, thereby simultaneously processing large-scale noise and sudden change noise without the limit of noise form.
The invention also provides a magnetotelluric denoising system, which comprises: the device comprises a magnetotelluric pulse signal measuring module, a sounding data segment obtaining module, an approximate entropy determining module of a sounding data segment, a sounding data segment determining module to be denoised, a matrix constructing module to be denoised, a singular value decomposition module, a field source component matrix obtaining module, a component data segment obtaining module, an approximate entropy determining module of a component data segment, a denoised sounding data segment obtaining module and a reconstructed denoised magnetotelluric pulse signal obtaining module.
And the magnetotelluric pulse signal measuring module is used for measuring magnetotelluric pulse signals of different depths of the earth by utilizing a magnetotelluric sounding method.
And the sounding data segment obtaining module is used for carrying out overlapping segmentation on the magnetotelluric pulse signals to obtain a plurality of sounding data segments.
And the approximate entropy determining module of the sounding data segment is used for determining the approximate entropy of each sounding data segment based on the approximate entropy theory.
And the depth data segment to be denoised determining module is used for determining the depth data segment with the approximate entropy smaller than or equal to the approximate entropy threshold as the depth data segment to be denoised and determining the depth data segment with the approximate entropy larger than the approximate entropy threshold as the effective depth data segment.
And the matrix construction module to be denoised is used for constructing the matrix to be denoised according to the sounding data segment to be denoised.
And the singular value decomposition module is used for performing singular value decomposition on the matrix to be denoised to obtain a singular value matrix, a left singular matrix and a right singular matrix of the matrix to be denoised.
And the field source component matrix obtaining module is used for multiplying each diagonal element of the singular value matrix, the column vector corresponding to the left singular matrix and the row vector corresponding to the right singular matrix to obtain a plurality of field source component matrixes of the sounding data section to be denoised.
A component data segment obtaining module for converting each field source component matrix into a component data segment.
And the approximate entropy determining module of the component data segments is used for determining the approximate entropy of each component data segment based on the approximate entropy theory.
And the denoised depth measurement data section obtaining module is used for screening the component data sections of which the approximate entropies are larger than that of the depth measurement data section to be denoised, and summing the screened component data sections to obtain the denoised depth measurement data section.
And the reconstruction denoising magnetotelluric pulse signal obtaining module is used for reconstructing the denoised depth sounding data segment and the effective depth sounding data segment to obtain a reconstruction denoising magnetotelluric pulse signal.
The sounding data segment obtaining module specifically comprises:
a sounding data segment obtaining submodule for obtaining magnetotelluric pulse signal x1,x2,…xwPerforming overlapping segmentation processing to obtain each depth measurement data segment Xv＝{x(v-1)*N+1,x(v-1)*N+2,…,x(v-1)*N+N}。
Wherein x is1,x2,…xwRespectively the 1 st, 2 nd and w th magnetotelluric pulse signals, w is the length of the magnetotelluric pulse signal, x(v-1)*N+1、x(v-1)*N+2、x(v-1)*N+NThe data are respectively the 1 st, 2 nd and nth data of the v-th sounding data segment, N is the sounding data segment length, v is 1,2, …, l is the number of sounding data segments, l is (w- α N)/N (1- α), and α is the data overlap ratio.
The approximate entropy determination module of the sounding data segment specifically comprises:
and the embedding dimension initialization submodule is used for initializing the embedding dimension h of the Hankel matrix.
And the Hankel matrix construction submodule is used for respectively constructing a Hankel matrix with an embedded dimension h for each depth sounding data segment. Each column in the hankel matrix is a mode vector.
Maximum distance determination submodule for utilizing a formula
And determining the maximum distance between the p-th mode vector in the Hankel matrix and the q-th mode vector of the Hankel matrix.
And the number counting submodule is used for counting the number of the p-th mode vector in the Hankel matrix, wherein the maximum distance between the p-th mode vector and all the mode vectors of the Hankel matrix is smaller than or equal to the mode fault-tolerant threshold.
An approximate probability value determining submodule for utilizing a formula C according to the number based on the approximate probability theoryp h(r)＝Np h(r)/(N-h +1), determining the approximate probability value of the p-th mode vector in the Hankel matrix.
An approximate probability average determination submodule for utilizing a formula based on the approximate probability value
And determining the approximate probability average value of each sounding data segment under the h embedding dimension.
And the circulation submodule is used for increasing the embedding dimension h by 1, returning to the step of constructing a Hankel matrix with the embedding dimension h for each sounding data segment respectively, and determining the approximate probability average value of each sounding data segment under the embedding dimension h + 1.
An approximate entropy determination submodule of the sounding data segment for utilizing a formula according to the approximate probability average value of each sounding data segment in the h embedding dimension and the approximate probability average value of each sounding data segment in the h +1 embedding dimension
An approximate entropy for each sounding data segment is determined.
Wherein d ismax[y(p),y(q)]Is the maximum distance between the p-th mode vector in the Heckel matrix and the q-th mode vector in the Heckel matrix, y (p) is the p-th mode vector in the Heckel matrix, y (q) is the q-th mode vector in the Heckel matrix, x (p + delta) is the delta-th element of the p-th mode vector in the Heckel matrix, x (q + delta) is the delta-th element of the q-th mode vector in the Heckel matrix, delta is 1,2, …, h, Cp h(r) is the approximate probability value of the p-th pattern vector in the Hankel matrix, Np h(r) is the number of maximum distances less than or equal to the mode fault tolerance threshold, r is the mode fault tolerance threshold, N is the sounding data segment length,
is an approximate probability average of the sounding data segment in the h embedding dimension,
is the approximate probability average, Ae, of the sounding data segment in the h +1 embedding dimensioniIs the approximate entropy of the ith sounding data segment.
The matrix construction module to be denoised specifically comprises:
an equal division module used for dividing the a-th depth sounding data segment X to be denoised according to the preset matrix dimensiona＝{x(a-1)*N+1,x(a-1)*N+2,…,x(a-1)*N+NAnd equally dividing into a plurality of sections.
A sub-module for constructing the matrix to be denoised according to each segment after being equally divided
Wherein, XaFor the a-th sounding data segment to be denoised, x(a-1)*N+1、x(a-1)*N+2、x(a-1)*N+n、x(a-1)*N+1+n、x(a-1)*N+2n、xa*N+1-nAnd x(a-1)*N+NRespectively representing the 1 st, 2 nd, nth, N +1 st, 2 nth, N-N +1 th and nth data of the depth measuring data segment to be denoised, wherein N is the length of the depth measuring data segment, N is the length of a preset matrix, N is N/m, and m is the dimension of the preset matrix.
The magnetotelluric denoising technology based on approximate entropy and SVD can reduce the phenomenon of 'over-processing' data under certain conditions; the requirement on the signal-to-noise ratio of data is low; the noise with complex forms can be suppressed.
The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other. For the system disclosed by the embodiment, the description is relatively simple because the system corresponds to the method disclosed by the embodiment, and the relevant points can be referred to the method part for description.
The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.
Claims (8)
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011309115.0A CN112285793B (en) | 2020-11-20 | 2020-11-20 | A kind of magnetotelluric denoising method and system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011309115.0A CN112285793B (en) | 2020-11-20 | 2020-11-20 | A kind of magnetotelluric denoising method and system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112285793A CN112285793A (en) | 2021-01-29 |
| CN112285793B true CN112285793B (en) | 2021-08-13 |
Family
ID=74398495
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011309115.0A Expired - Fee Related CN112285793B (en) | 2020-11-20 | 2020-11-20 | A kind of magnetotelluric denoising method and system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112285793B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113568058B (en) * | 2021-07-20 | 2022-06-03 | 湖南师范大学 | A method and system for separation of magnetotelluric signal-to-noise based on multi-resolution singular value decomposition |
| CN114239651B (en) * | 2021-12-15 | 2024-10-18 | 吉林大学 | Earth electromagnetic denoising method and system |
| CN119758800B (en) * | 2024-11-19 | 2025-07-15 | 江苏赫玛信息科技有限公司 | Intelligent manufacturing process data acquisition method and system |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105467460A (en) * | 2015-12-28 | 2016-04-06 | 中国石油天然气集团公司 | Method and device for electromagnetic prospecting |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090265111A1 (en) * | 2008-04-16 | 2009-10-22 | Kjt Enterprises, Inc. | Signal processing method for marine electromagnetic signals |
| CN103584872B (en) * | 2013-10-29 | 2015-03-25 | 燕山大学 | Psychological stress assessment method based on multi-physiological parameter fusion |
| CN103837898B (en) * | 2014-02-24 | 2016-08-17 | 吉林大学 | High-density electric near-end dipole electromagnetic sounding method |
| CN111709279B (en) * | 2020-04-30 | 2023-07-21 | 天津城建大学 | An Algorithm for Separating Microseismic Noise Mixed Signals Using SVD-EMD Algorithm |
-
2020
- 2020-11-20 CN CN202011309115.0A patent/CN112285793B/en not_active Expired - Fee Related
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105467460A (en) * | 2015-12-28 | 2016-04-06 | 中国石油天然气集团公司 | Method and device for electromagnetic prospecting |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112285793A (en) | 2021-01-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Dong et al. | Non-iterative denoising algorithm for mechanical vibration signal using spectral graph wavelet transform and detrended fluctuation analysis | |
| CN110865357B (en) | Laser radar echo signal noise reduction method based on parameter optimization VMD | |
| CN112285793B (en) | A kind of magnetotelluric denoising method and system | |
| Saad et al. | Unsupervised deep learning for single-channel earthquake data denoising and its applications in event detection and fully automatic location | |
| CN108345033B (en) | A first-arrival detection method for microseismic signals in time-frequency domain | |
| CN113568058B (en) | A method and system for separation of magnetotelluric signal-to-noise based on multi-resolution singular value decomposition | |
| CN110646851B (en) | An adaptive threshold seismic random noise suppression method based on Shearlet transform | |
| CN110031899A (en) | Compressed sensing based weak signal extraction algorithm | |
| CN105676292A (en) | 3D earthquake data de-noising method based on 2D curvelet transform | |
| CN107179550B (en) | A kind of seismic signal zero phase deconvolution method of data-driven | |
| CN109633761B (en) | A method for reducing power frequency noise of magnetic resonance signal based on wavelet transform modulus maximum method | |
| CN108961181A (en) | A kind of ground penetrating radar image denoising method based on shearlet transformation | |
| Il et al. | An appropriate thresholding method of wavelet denoising for dropping ambient noise | |
| CN114091538B (en) | A discriminant loss convolutional neural network intelligent denoising method based on signal characteristics | |
| CN106405504B (en) | The Coherent Noise in GPR Record denoising method of combined shearing wave conversion and singular value decomposition | |
| CN117076858A (en) | A method and system for suppressing strong low-frequency electromagnetic interference based on deep learning | |
| CN107015124B (en) | A Method of Partial Discharge Signal Interference Suppression Based on Frame Adaptive Sparse Decomposition | |
| CN116520317B (en) | Ground Penetrating Radar Signal Denoising Method Combining Two-Dimensional VMD and DT-CWT | |
| Wang et al. | A joint framework for seismic signal denoising using total generalized variation and shearlet transform | |
| CN112817056A (en) | Magnetotelluric signal denoising method and system | |
| CN116401509A (en) | Underwater ultrasonic signal denoising method based on improved EWT-FastICA combination | |
| CN113655534A (en) | Noise suppression method of NMR FID signal based on multilinear singular value tensor decomposition | |
| Kim et al. | Statistical modeling and denoising of microseismic signal for dropping ambient noise in wavelet domain | |
| CN109782346B (en) | A method of collecting footprints based on morphological component analysis | |
| CN104765069A (en) | Method of suppressing adjacent shot interference in case of synchronous excitation acquisition |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210813 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |