CN109101954B - Method for removing tidal strain components of borehole strain data based on minimum noise separation - Google Patents

Abstract

The invention relates to a method for removing tidal strain components of borehole strain data based on minimum noise separation, which adopts a self-adaptive window width filtering algorithm to estimate high-frequency signals of a borehole strain data sample matrix; calculating a covariance matrix of the estimated high-frequency signal, and decomposing an eigenvalue; constructing an adjustment matrix by using the covariance matrix of the high-frequency signals; adjusting the covariance matrix of the drilling strain data sample matrix by using the constructed adjustment matrix; carrying out minimum noise separation transformation on the covariance matrix of the adjusted borehole strain data sample matrix to obtain a minimum noise separation transformation component, wherein the matrix formed by the minimum noise separation transformation component is represented as; calculating a harmonic period of the minimum noise separation component; selecting a minimum noise separation component corresponding to the harmonic period of the tidal strain component for reconstruction; and calculating the data for removing the tidal strain component, and specifically removing the tidal strain component without calculating a theoretical solid tide.

Description

Method for removing tidal strain components of borehole strain data based on minimum noise separation

Technical Field

The invention relates to the field of borehole strain data processing, in particular to a method for removing tidal strain components of borehole strain data based on minimum noise separation.

Background

The deformation of the earth surface, the movement of the structure in the earth crust and various geological disasters generated by the movement are closely related to the action of the stress of the earth crust, and the change of the stress-strain state of the earth crust is the most direct cause of fracture, wrinkle and even earthquake. Borehole strain observation is an important earthquake precursor observation means by accurately observing the continuous change of the internal strain state of the stratum according to time, and finding and mastering the time-space distribution and development change rules of the (long) medium term-short term-impending earthquake and post-earthquake adjustment of earthquake strain precursors. The solid tidal response can be clearly recorded due to the high precision of borehole strain observation. The solid tidal response will cause the borehole strain data to produce a tidal strain component that masks the short-term high frequency information of the crustal strain that causes the formation earthquake. The minimum noise separation method can theoretically judge the intrinsic dimensionality of data, separate signals in the data and solve the problem of removing tidal strain components in the borehole strain data.

CN106918836A discloses a borehole strain data anomaly extraction method based on principal component analysis, which removes periodic responses of borehole strain data by using a harmonic analysis method, and respectively represents weak changes of the crust by using eigenvalues and eigenvectors in the principal component analysis.

CN104390733A discloses a determination method of magnitude and direction of ground stress, wherein three strain patterns are arranged on an epoxy resin hollow bag body for measuring stress along the circumferential direction, the included angle of adjacent strain patterns is 120 degrees, each strain pattern is provided with four strain sheets, and the strain sheets on each strain pattern are sequentially rotated by 45 degrees; by utilizing 12 strain values in different directions in the circumferential direction, the ground stress and the stress direction can be accurately and effectively calculated.

The method comprises the steps of extracting high-frequency components of borehole strain data by a high-pass filtering method, and the like (extracting strain abnormality of Wenchuan earthquake body of Ning shan platform by wavelet-transfinite rate analysis, 2012). And analyzing the Wenchuan earthquake stress anomaly of the Ningshan platform by using a wavelet-transfinite rate method. Li cheng wu et al (preliminary analysis of surface strain tidal factor observed by borehole strain gauges, 2014) calculates borehole strain differential data, and analyzes the tidal factor change characteristics and the factors affecting the tidal factor dynamic change. However, no report has been made to remove the tidal strain component of borehole strain data by using a minimum noise separation method. The existing method does not consider the physical significance of each harmonic wave and lacks pertinence; some need to calculate theoretical solid tide, and the calculation amount is large.

Disclosure of Invention

The invention aims to provide a method for removing tidal strain components of borehole strain data based on minimum noise separation, which is used for removing the tidal strain components in the borehole strain data in a targeted manner.

The invention is realized by a method for removing tidal strain components of borehole strain data based on minimum noise separation, which comprises the following steps:

a. making the borehole strain data into a sample matrix Y (n x 1440);

b. estimating a high-frequency signal H (n x 1440) by adopting an adaptive window width filtering algorithm on a borehole strain data sample matrix Y (n x 1440);

c. calculating covariance matrix C of estimated high-frequency signal H (n x 1440)_H(n × n), and performing eigenvalue decomposition;

d. covariance matrix C using high frequency signals_H(n × n) to construct an adjustment matrix P;

e. adjusting the covariance matrix CY (n multiplied by n) of the borehole strain data sample matrix by using the constructed adjustment matrix P to obtain the covariance matrix C of the adjusted borehole strain data sample matrix_P(n×n)；

f. Covariance matrix C of adjusted borehole strain data sample matrix_P(n × n) performing minimum noise separation transformation to obtain a minimum noise separation transformation component, wherein a matrix formed by the minimum noise separation transformation component is represented by ψ;

g. calculating a harmonic period ψ (ω) of the minimum noise separation component;

h. selecting a minimum noise separation component corresponding to the harmonic period of the tidal strain component for reconstruction;

i. data for removing the tidal strain component is calculated, and the result is output.

Further, the step a comprises: selecting a group of drilling strain minute data, and expressing the drilling strain data as follows according to a time sequence:

X₁＝[X₁(1)，X₁(2)，...，X₁(1440)]，...，X_n＝[X_n(1)，X_n(2)，...，X_n(1440)]obtaining a sample matrix y ═ X₁，X₂，X₃，...X_n]^TThe expression of the sample matrix Y (n × 1440) is:

further, step b includes filtering the prepared sample matrix Y (n × 1440) by an adaptive window width filter, and estimating a high-frequency signal in the sample matrix, where the adaptive window width filter is expressed as follows:

wherein, O is filtered data, and Y is input sample data; w is the adaptive window width, and the expression is as follows:

wherein, W_LIs a selectable minimum window width, W_UΔ (j) is the second order difference of the sample matrix y (n x 1440) for the selectable maximum window width;

the estimated high frequency signal is then:

H＝Y-O (4)

where H denotes H (n × 1440) as the estimated high frequency signal, Y denotes Y (n × 1440) as the borehole strain data sample matrix, and O denotes O (n × 1440) as the adaptive window width filtered data.

Further, step C includes calculating a covariance matrix C of the estimated high-frequency signal H (n × 1440)_H(n × n), covariance matrix C_H(n × n) element γ_pqCalculated by the formula (5),

wherein the content of the first and second substances,

and

row i, minutes p and minutes q, respectively;

and

the mean values of the p-th and q-th minute data of the N rows of data respectively;

for the calculated covariance matrix C_H(n × n) characteristic value decomposition is performed, and the expression is as follows:

in the formula, D_HA diagonal matrix with eigenvalues arranged in descending order; r_HIs the eigenvector matrix corresponding to the eigenvector matrix, and T is the transpose of the matrix.

Further, step d includes utilizing the covariance matrix C of the estimated high frequency signal H (n x 1440)_H(n × n) constructing an adjustment matrix P such that the adjustment matrix satisfies P^TC_HWhere E is an identity matrix, the expression of the adjustment matrix is as follows:

P＝R_HD_H ^-0.5 (7)

wherein D is_HA diagonal matrix with eigenvalues arranged in descending order; r_HAnd the characteristic vector matrix is corresponding to the characteristic value diagonal matrix.

Further, step e includes utilizing the constructed adjustment matrix P to a covariance matrix C of the borehole strain data sample matrix_Y(n × n) is adjusted by the expression: c_P＝P^TC_YP, wherein C_P(n) is a covariance matrix of the adjusted borehole strain data sample matrix.

Further, step f includes a covariance matrix C of the adjusted borehole strain data sample matrix_P(n × n) is subjected to eigenvalue decomposition, and the expression is as follows:

wherein D is_PIs a diagonal matrix with eigenvalues of covariance matrix of adjusted borehole strain data sample matrix arranged in descending order, V_PThe characteristic vector matrix corresponding to the characteristic value of the covariance matrix of the adjusted borehole strain data sample matrix is obtained, and the transformation matrix R of minimum noise separation is expressed as: r ═ PV_PWhere P is the adjustment matrix, the matrix ψ formed by the minimum noise separation transformation into bins can be expressed as:

ψ＝YR＝Y(PV_P)＝[ψ₁，ψ₂，...，ψ_n]，n＝1，2，...，m， (9)

where Y denotes Y (n x 1440) as a borehole strain data sample matrix, #₁，ψ₂，...，ψ_nIs the minimum noise separation component.

Further, step g includes using fourier transform to find the harmonic period of each minimum noise separation component, which is calculated as follows:

wherein psi_n(t) is a minimum noise separation component, and ψ (ω) is a harmonic period of the minimum noise separation component.

Further, the step h comprises selecting k minimum noise separation components corresponding to the harmonic period of the tidal strain component and the harmonic period of the harmonic period, and reconstructing, wherein the reconstructed data expression is as follows:

Y_MNF＝ψR^-1＝[ψ₁，ψ₂，...，ψ_k，0，0，...，0]R^-1，1≤k≤m， (11)。

further, step i includes transforming the reconstructed data matrix Y (n x 1440) with the borehole strain data sample matrix Y and the minimum noise separation transform_MNF(n x 1440) to remove tidal StrainData Y of the component_S(n × 1440), the expression of which is as follows:

Y_S＝Y-Y_MNF (12)

where Y (n x 1440) is the borehole strain data sample matrix, Y_MNF(n × 1440) is a data matrix reconstructed by the minimum noise separation transform.

Compared with the prior art, the invention has the beneficial effects that:

the method can determine the intrinsic dimension of the data, and remove the tidal strain component in a targeted manner by utilizing the characteristic of the tidal strain component in the borehole strain data on the frequency without calculating the theoretical solid tide.

Drawings

FIG. 1 is a flow chart of a method for removing tidal strain components of borehole strain data based on minimum noise separation;

FIG. 2 is a sample matrix diagram;

FIG. 3 is a minimum noise separation component diagram;

FIG. 4 (a), (b), (c), (d), (e), (f), (g), (h), (i) are spectrograms of the 9 least noise separated components corresponding to the harmonic period of the tidal strain component;

FIG. 5 is a graph of data reconstructed using selected minimum noise separation components;

FIG. 6 is a plot of time series data with the tidal strain component removed.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

Take the drilling strain minute data of Guzan earthquake precursor monitoring station in Sichuan as an example. The data were measured using an YRY quarter-gauge borehole strain gauge from 1/2012 to 12/31/2013.

a, recording the drilling strain data S₁The daily data for borehole strain is expressed as: x₁＝[X₁(1)，X₁(2)，...，X₁(1440)]，...，X₇₃₁＝[X₇₃₁(1)，X₇₃₁(2)，...，X₇₃₁(1440)]A sample matrix y ═ X can be obtained₁，X₂，...，X₇₃₁]^TAs shown in fig. 2. The expression of the sample matrix y (731 × 1440) is:

and b, performing high-frequency signal estimation on the borehole strain data sample matrix by adopting an adaptive window width filtering algorithm, namely filtering the prepared sample matrix Y (731 multiplied by 1440) by using an adaptive window width filter to estimate high-frequency signals in the sample matrix. The expression for the adaptive window width filter is as follows:

where O is filtered data and Y is input sample data. W is the adaptive window width, and the expression is as follows:

wherein, W_LIs a selectable minimum window width, W_UIs the maximum window width that can be selected. Δ (j) is the second order difference of the sample matrix Y (731 × 1440).

The estimated high frequency signal is then:

H＝Y-O (4)

where H (731X 1440) is the estimated high frequency signal, Y (731X 1440) is the borehole strain data sample matrix, and O (731X 1440) is the adaptive window width filtered data.

C, calculating covariance matrix C of estimated high-frequency signal H (731X 1440)_H(731X 731), covariance matrix C_H(731X 731) element gamma_pqCalculated by the formula (5),

wherein the content of the first and second substances,

and

row i, minutes p and minutes q, respectively;

and

are the mean values of the p and q minute data of the N line data, respectively.

For the calculated covariance matrix C_H(731 × 731) performing eigenvalue decomposition, the expression is as follows:

in the formula, D_HA diagonal matrix with eigenvalues arranged in descending order; r_HIs the eigenvector matrix corresponding to the eigenvector matrix, and T is the transpose of the matrix.

d, constructing an adjustment matrix, and using the covariance matrix C of the estimated high-frequency signal H (731X 1440)_H(731 × 731) constructing an adjustment matrix P such that the adjustment matrix satisfies P^TC_HP ═ E, where E is the identity matrix. The expression of the adjustment matrix is as follows:

P＝R_HD_H ^-0.5 (7)

wherein D is_HA diagonal matrix with eigenvalues arranged in descending order; r_HAnd the characteristic vector matrix is corresponding to the characteristic value diagonal matrix.

e, using the covariance matrix C of the constructed adjustment matrix P to the borehole strain data sample matrix_Y(731 × 731) adjusting, by the expression: c_P＝P^TC_YP, wherein C_P(731X 731) is the adjusted borehole strain data sample covariance matrix, C_YThe (731 × 731) can be calculated by the formula (5).

f, covariance matrix C of adjusted borehole strain data sample matrix_P(731 × 731) performing eigenvalue decomposition, wherein the expression is as follows:

wherein D is_PIs a diagonal matrix with eigenvalues of covariance matrix of adjusted borehole strain data sample matrix arranged in descending order, V_PAnd the characteristic vector matrix is corresponding to the characteristic value of the covariance matrix of the adjusted drilling strain data sample matrix. The transformation matrix R for minimum noise separation can be expressed as: r ═ PV_PWhere P is the adjustment matrix. The matrix ψ formed by the minimum noise separation transformation into bins can be expressed as:

ψ＝YR＝Y(PV_P)＝[ψ₁，ψ₂，...，ψ_n]，n＝1，2，...，731， (9)

where Y (731X 1440) is the borehole strain data sample matrix, #₁，ψ₂，...，ψ_nIs the minimum noise separation component. Fig. 3 is a minimum noise separation component diagram.

g, calculating the harmonic period of each minimum noise separation component by utilizing Fourier transform, wherein the calculation formula is as follows:

wherein psi_n(t) is a minimum noise separation component, and ψ (ω) is a harmonic period of the minimum noise separation component.

h, selecting 9 minimum noise separation components corresponding to the harmonic period of the tidal strain component to reconstruct, wherein FIG. 4 is a frequency spectrum diagram of the 9 minimum noise separation components corresponding to the harmonic period of the tidal strain component. The reconstructed data expression is as follows:

Y_MNF＝ψR^-1＝[ψ₁，ψ₂，...，ψ₉，0，0，...，0]R^-1， (11)

the reconstructed data is shown in fig. 5.

i, finding data Y from which the tidal Strain component has been removed_S(731 × 1440), the expression of which is as follows:

Y_S＝Y-Y_MNF (12)

where Y (731X 1440) is the borehole strain data sample matrix, Y_MNF(731 × 1440) is the data matrix reconstructed by the minimum noise separation transform. Y is_S(731 × 1440) is data from which the tidal strain component was removed. The time series data with the tidal strain component removed is shown in fig. 6.

The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and principle of the present invention are intended to be included within the scope of the present invention.

Claims (7)

1. A method for removing tidal strain components of borehole strain data based on minimum noise separation, the method comprising:

a. making the borehole strain data into a sample matrix Y (n x 1440);

b. estimating a high-frequency signal H (n x 1440) by adopting an adaptive window width filtering algorithm on a borehole strain data sample matrix Y (n x 1440);

c. calculating covariance matrix C of estimated high-frequency signal H (n x 1440)_H(n × n), and performing eigenvalue decomposition;

d. covariance matrix C using high frequency signals_H(n × n) to construct an adjustment matrix P;

e. covariance matrix C of borehole strain data sample matrix using constructed adjustment matrix P_Y(n x n) to obtain covariance matrix C of adjusted borehole strain data sample matrix_P(n×n)；

f. Covariance matrix C of adjusted borehole strain data sample matrix_P(n × n) performing minimum noise separation transformation to obtain a minimum noise separation transformation component, wherein a matrix formed by the minimum noise separation transformation component is represented by ψ;

g. calculating a harmonic period ψ (ω) of the minimum noise separation component;

h. selecting a minimum noise separation component corresponding to the harmonic period of the tidal strain component for reconstruction;

i. calculating data for removing the tidal strain component, and outputting a result;

step g comprises the following steps of utilizing Fourier transform to solve the harmonic period of each minimum noise separation component, wherein the calculation formula is as follows:

wherein psi_n(t) is a minimum noise separation component, and ψ (ω) is a harmonic period of the minimum noise separation component;

step h comprises selecting k minimum noise separation components corresponding to the harmonic period of the tidal strain component and the harmonic period of the tidal strain component for reconstruction, wherein the data expression after reconstruction is as follows:

Y_MNF＝ψR^-1＝[ψ₁，ψ₂，...，ψ_k，0，0，...，0]R^-1，1≤k≤m，

step i includes transforming the reconstructed data matrix Y (n x 1440) with the borehole strain data sample matrix Y and minimum noise separation_MNF(n x 1440) data Y for removing tidal Strain component_S(n × 1440), the expression of which is as follows:

Y_S＝Y-Y_MNF

where Y (n x 1440) is the borehole strain data sample matrix, Y_MNF(n × 1440) is a data matrix reconstructed by the minimum noise separation transform.

2. The method of claim 1, wherein step a comprises: selecting a group of drilling strain minute data, and expressing the drilling strain data as follows according to a time sequence:

X₁＝[X₁(1)，X₁(2)，...，X₁(1440)]，…，X_n＝[X_n(1)，X_n(2)，...，X_n(1440)]obtaining a sample matrix Y ═ X₁，X₂，X₃，...X_n]^TThe expression of the sample matrix Y (n × 1440) is:

3. the method of claim 1, wherein step b comprises filtering the prepared sample matrix Y (nx 1440) through an adaptive window width filter to estimate the high frequency signal in the sample matrix, the adaptive window width filter being expressed as follows:

wherein, O is filtered data, and Y is input sample data; w is the adaptive window width, and the expression is as follows:

wherein, W_LIs a selectable minimum window width, W_UΔ (j) is the second order difference of the sample matrix Y (n × 1440) for the selectable maximum window width;

the estimated high frequency signal is then:

H＝Y-O (4)

where H denotes H (n × 1440) as the estimated high frequency signal, Y denotes Y (n × 1440) as the borehole strain data sample matrix, and O denotes O (n × 1440) as the adaptive window width filtered data.

4. The method of claim 1, wherein step C comprises calculating a covariance matrix C of the estimated high-frequency signal H (nx 1440)_H(n × n), covariance matrix C_H(n × n) element γ_pqCalculated by the formula (5),

wherein the content of the first and second substances,

and

row i, minutes p and minutes q, respectively;

and

the mean values of the p-th and q-th minute data of the N rows of data respectively;

for the calculated covariance matrix C_H(n × n) characteristic value decomposition is performed, and the expression is as follows:

in the formula, D_HA diagonal matrix with eigenvalues arranged in descending order; r_HIs an eigenvector matrix corresponding to the eigenvector matrix,

is R_HTransposing of the matrix.

5. The method of claim 1, characterized in thatCharacterized in that step d comprises using the covariance matrix C of the estimated high-frequency signal H (n x 1440)_H(n × n) constructing an adjustment matrix P such that the adjustment matrix satisfies P^TC_HWhere E is an identity matrix, the expression of the adjustment matrix is as follows:

P＝R_HD_H ^-0.5 (7)

wherein D is_HA diagonal matrix with eigenvalues arranged in descending order; r_HAnd the characteristic vector matrix is corresponding to the characteristic value diagonal matrix.

6. The method of claim 1, wherein step e comprises using the covariance matrix C of the constructed adjustment matrix P versus the borehole strain data sample matrix_Y(n × n) is adjusted by the expression: c_P＝p^TC_YP, wherein C_P(n) is a covariance matrix of the adjusted borehole strain data sample matrix.

7. The method of claim 1, wherein step f comprises a covariance matrix C of the adjusted borehole strain data sample matrix_P(n × n) is subjected to eigenvalue decomposition, and the expression is as follows:

wherein D is_PIs a diagonal matrix with eigenvalues of covariance matrix of adjusted borehole strain data sample matrix arranged in descending order, V_PThe characteristic vector matrix corresponding to the characteristic value of the covariance matrix of the adjusted borehole strain data sample matrix is obtained, and the transformation matrix R of minimum noise separation is expressed as: r ═ PV_PWhere P is the adjustment matrix, the matrix ψ formed by the minimum noise separation transformation into bins can be expressed as:

ψ＝YR＝Y(PV_P)＝[ψ₁，ψ₂，...，ψ_n]，n＝1，2，...，m， (9)

wherein, Y tableY (n x 1440) is shown as the borehole strain data sample matrix, #₁，ψ₂，...，ψ_nIs the minimum noise separation component.

Images