A millennium-scale period identification method and device based on bispectrum analysis

Abstract

The present invention discloses a method and device for identifying millennium-scale cycles based on bispectral analysis, including: (1) selecting paleoclimate proxy indicators, performing power spectrum analysis after preprocessing, and identifying thickness cycles corresponding to astronomical orbital cycles; (2) filtering the frequencies corresponding to the thickness cycles, selecting the filtered curve to establish a time-depth conversion model to obtain the time-domain sequence of the proxy indicator; (3) performing power spectrum analysis on the time-domain sequence to identify astronomical orbital cycles and millennium-scale cycles; (4) importing the time-domain sequence into the Matlab high-order spectral analysis toolbox, and calculating the non-linear phase coupling relationship between two sets of frequencies through bispectral analysis; (5) exporting the bispectral diagram of the entire time-domain curve in the time domain in step (4), observing and identifying strong response regions with high correlation, and further verifying the existence of the millennium-scale cycle. The present invention can solve the problem of identifying millennium-scale cycles and short-scale climate change cycles.

Description

Technical Field

[0001] The present invention belongs to the technical field of cyclostratigraphic information processing, and in particular relates to a method and device for identifying millennial-scale cycles based on bispectral analysis. Background Art

[0002] In recent years, with the improvement of XRF element scanning technology, XRF element scanning of cores has become a means to easily obtain a large amount of high-resolution geochemical data, which provides an opportunity for geologists and astronomers to detect astronomical orbital cycles preserved in sedimentary rocks.

[0003] For example, the Chinese patent document with the publication number CN104089964A discloses a dating method for logging Milankovitch cycle analysis, including the following steps: (1) obtaining sedimentary record data and corresponding depth values of sediments at different depths in the target interval; (2) obtaining core samples and depth values of sediments at different depths in the target interval, and analyzing the core samples to determine the depths where event sediments are located; (3) removing the sedimentary record data at the depths where event sediments are located from the sedimentary record data to obtain data and depth values to be subjected to spectral analysis; (4) matching the data and depth values to be subjected to spectral analysis with pre-obtained astronomical cycles to obtain the time when the sediments were formed.

[0004] The Chinese patent document with the publication number CN114578444A discloses an astronomical stratigraphic cycle division method based on wavelet multi-scale analysis, including the following steps: Step 1. Comprehensively compare the characteristics of various logging curves and select sensitive logging curves; Step 2. Perform spectral analysis on the logging data to find the dominant frequencies controlled by orbital parameters; Step 3. According to the time control points of known formation ages, obtain the theoretical curves of astronomical orbits in the geological history period of the target interval; Step 4. Extract the dominant scale factors through power spectra and perform multi-scale wavelet decomposition on the logging curves at this scale value; Step 5. Compare the wavelet decomposition results with the theoretical curves of astronomical orbits and select the target curves that can be used for tuning to divide astronomical stratigraphic cycles.

[0005] However, the hundred-thousand-year-scale cycles (astronomical orbital cycles) interfere with the changes of short cycles to a certain extent, making it difficult to identify millennial-scale cycles. Therefore, there is an urgent need to design a method for identifying millennial-scale cycles. Summary of the Invention

[0006] The present invention provides a method for identifying millennial-scale cycles based on bispectral analysis, which can solve the problems of identifying millennial-scale cycles and short-scale climate change cycles.

[0007] A method for identifying millennial-scale cycles based on bispectral analysis, characterized by including:

[0008] (1) Select paleoclimate proxy indicators for detecting millennial-scale periodic climate changes. After preprocessing (such as interpolation and removing background white noise), perform power spectrum analysis to identify thickness cycles corresponding to astronomical orbital periods.

[0009] (2) Perform Gaussian band-pass filtering on the frequencies corresponding to the identified thickness cycles, select the filtered curve to establish a time-depth conversion model, and obtain the time-domain sequence of the proxy indicator based on the time-depth conversion model.

[0010] (3) Perform power spectrum analysis on the time-domain sequence to identify astronomical orbital periods and millennial-scale periods in the time-domain sequence.

[0011] (4) Import the time-domain sequence into the Matlab high-order spectral analysis toolbox. Through bispectral analysis, couple the phases of the astronomical orbital periods in the time-domain sequence at frequencies f1, f2, and f3 respectively. The energy is non-linearly transferred between these frequencies and redistributed on the frequency spectrum, forming a skewed and / or asymmetric cyclic geometric structure.

[0012] Among them, frequencies f1 and f2 correspond to the x-axis and y-axis respectively, and f3 = f1 + f2.

[0013] (5) Compare and identify the skewed and / or asymmetric geometric structures formed in the bispectral analysis, export the bispectral diagram of the entire time-domain curve in the time domain in step (4), and use the x-axis, y-axis, and projection diagonal line to observe and identify the strongly responsive regions with high correlation. Combine with the millennial-scale period identified in step (3) to further verify the existence of this millennial-scale period.

[0014] In the present invention, when the phases of the bispectrum are coupled, the energy is non-linearly transferred between these frequencies and redistributed on the frequency spectrum; this energy transfer forms lower and higher harmonics and forms a skewed and (or) asymmetric cyclic geometric structure, and stronger periodicity can be found in the bispectrum of the selected time interval. Traditional spectral analysis methods ignore the astronomical cycle phase information hidden in the signal, while the method of the present invention detects second-order phase coupling with reasonable time resolution, which is helpful for the identification of millennial-scale periods.

[0015] In step (1), the paleoclimate proxy indicators for detecting millennial-scale periodic climate changes include major elements, trace elements, or their ratios obtained by XRF element scanning.

[0016] Astronomical orbital periods refer to the long eccentricity period (405 thousand years), short eccentricity period (about 100 thousand years), slope period (about 40 thousand years), and precession period (about 20 thousand years). Millennial-scale periods refer to periods less than the precession period and greater than 1 thousand years.

[0017] Based on the ratio relationship of eccentricity, slope, and precession period, the corresponding thickness periods are identified; the relationship between the thickness period and the corresponding depth should satisfy the estimation of the sedimentation rate of the regional geological background.

[0018] In step (2), the frequency corresponding to the thickness period is the reciprocal of the thickness.

[0019] In step (3), identifying the astronomical orbital periods in the time domain sequence specifically includes: the identified frequencies corresponding to the long eccentricity, short eccentricity, slope, and precession period are 0.00245, 0.001, 0.0025, and 0.005 cycles / kiloyear respectively.

[0020] In step (4), assuming that the autocorrelation function of the time domain sequence {x(n)} is expressed as:

[0021] r(τ) = E{x(n)x(n + τ)}

[0022] The third-order matrix formula of {x(n)} is as follows:

[0023] R(τ1, τ2) = E{x(n)x(n + τ1)x(n + τ2)}

[0024] The bispectrum is expressed in the form of a frequency characteristic function and is the two-dimensional Fourier transform of R(τ1, τ2). The formula is as follows:

[0025] B(f1, f2) = E{X(f1)X(f2)X * (f1 + f2)}

[0026] In the formula, τ1, τ2 are random times, E{} is the mathematical expectation, X(f) is the Fourier transform of {x(n)}, X*(f1 + f2) is the complex conjugate; B(f1, f2) reflects the non-linear interaction of the astronomical signal components in the data at f1 and f2, and is considered to be the decomposition of the skewness of the signal in the frequency domain range;

[0027] The bispectral coherence is used to describe the asymmetric phase coupling of the signal. The above bispectrum is normalized to obtain the bispectral coherence, as shown in the formula:

[0028]

[0029] In the formula, P(ω1), P(ω2), P(ω1 + ω2) are the values of the power spectrum of x(n) at ω1, ω2, and ω1 + ω2.

[0030] In step (5), if there is no millennium-scale periodic signal at this frequency, the bispectrum coefficient is close to 0. When a millennium-scale periodic signal is detected, the bispectrum coefficient no longer remains zero. The more significant the millennium-scale periodic signal is, the closer the bispectrum coefficient is to 1. Further, the frequency corresponding to this period is determined according to the spectral peak on the bispectrum diagram, and the reciprocal of the frequency is the millennium-scale period.

[0031] A millennium-scale period recognition device based on bispectrum analysis includes a memory and one or more processors. An executable code is stored in the memory. When the one or more processors execute the executable code, it is used to implement the above-mentioned millennium-scale period recognition method.

[0032] Compared with the prior art, the present invention has the following beneficial effects:

[0033] The traditional Fourier transform (power spectrum analysis) is a method for analyzing stationary signals. However, most of the astronomical periodic signals recorded in geology have the characteristics of non-stationarity, and their frequency, amplitude, and phase will change with time. Bispectrum analysis is a method that can analyze phase information and can detect the quadratic phase coupling in signals. For high-resolution climate proxy indicators, power spectrum analysis cannot distinguish between millennium-scale periodic signals and noise signals, and the extracted frequency may not belong to the millennium-scale climate change cycle. However, the present invention uses bispectrum analysis to identify the associations and couplings between the frequency components of known astronomical orbital parameter periods (long eccentricity, short eccentricity, slope, and precession), so as to identify the millennium-scale climate change cycle. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a schematic flow chart of a millennium-scale period recognition method based on bispectrum analysis provided by the present invention;

[0035] Figure 2 It is a spectrum analysis diagram of the geochemical element data curve of a certain drilling well after pretreatment in an embodiment of the present invention;

[0036] Figure 3 It is the Gaussian band-pass filtering process in the depth domain of a certain drilling well in an embodiment of the present invention;

[0037] Figure 4 It is a time-depth conversion model constructed with a 100,000-year corresponding period for a certain drilling well in an embodiment of the present invention;

[0038] Figure 5 It is a bispectrum analysis diagram of a certain drilling well in an embodiment of the present invention, forming a skewed and / or asymmetric cyclic geometric structure map. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0039] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0040] It should be noted that, without conflict, the features in the following embodiments and implementation manners can be combined with each other.

[0041] As Figure 1 shown, a method for identifying millennium-scale cycles based on bispectral analysis includes the following steps:

[0042] S101, perform equidistant and background white noise removal processing on the selected geochemical element logging data, import it into open-source software for spectral analysis, and identify the corresponding astronomical orbital periods. As Figure 2 shown, the numbers in the figure are the reciprocals of the corresponding abscissas (frequencies), that is, the thickness periods.

[0043] According to the ratio relationship of each period of the astronomical orbital parameters and the ratio relationship of each thickness, it can be inferred that the 0.294m thickness period corresponds to the eccentricity period (100 thousand years).

[0044] S102, perform Gaussian band-pass filtering on the depth corresponding to the astronomical orbital period identified, that is, 0.294m. The center frequency of the filtering is selected as 1 / 0.294, that is, 3.4m -1 ( Figure 3 )

[0045] Establish a time-depth conversion model ( Figure 4 ) according to this filtering curve, that is, one depth period corresponds to 100 thousand years. Based on this time-depth conversion model, obtain the time-domain sequence of the geochemical element curve.

[0046] S103, perform power spectrum analysis (fast Fourier transform) on the time-domain index sequence, and identify the astronomical orbital periods (eccentricity, slope, and precession) and millennium-scale periods (less than half precession period, greater than one thousand-year period) in the time-domain sequence.

[0047] S104. Import the processed data into the Matlab High-Order Spectral Analysis (HOSA) toolbox. Through bispectral analysis, place the astronomical cycles in the time and index data sequences on the x-axis (frequency f1) and y-axis (frequency f2) and f3 (f3 = f1 + f2) respectively for phase coupling. Using the X-axis and Y-axis as references, form an inclined and / or asymmetric periodic geometric feature composed of the periodic coupling reflection intensity. The magnitude of energy transfer will be represented by the color intensity in the image and form the geometric image of each period in the coordinate axes.

[0048] Suppose the autocorrelation function of the time-domain sequence {x(n)} is expressed as:

[0049] r(τ) = E{x(n)x(n + τ)}

[0050] The third-order matrix formula of {x(n)} is as follows:

[0051] R(τ1, τ2) = E{x(n)x(n + τ1)x(n + τ2)}

[0052] The bispectrum is expressed in the form of a frequency characteristic function and is the two-dimensional Fourier transform of R(τ1, τ2). The formula is as follows:

[0053] B(f1, f2) = E{X(f1)X(f2)X * (f1 + f2)}

[0054] In the formula, τ1, τ2 are random times, E{} is the mathematical expectation, X(f) is the Fourier transform of {x(n)}, and X * (f1 + f2) is the complex conjugate; B(f1, f2) reflects the nonlinear interaction of the astronomical signal components in the data at f1 and f2 and is considered the decomposition of the skewness of the signal in the frequency domain range;

[0055] The bispectral coherence is used to describe the asymmetric phase coupling of the signal. Normalize the above bispectrum to obtain the bispectral coherence, as shown in the formula:

[0056]

[0057] In the formula, P(ω1), P(ω2), P(ω1 + ω2) are the values of the power spectrum of x(n) at ω1, ω2, and ω1 + ω2.

[0058] Such as Figure 5As shown, the color of the bispectrum shows the direction of energy transfer. White indicates that the spectral power transfers from two frequencies f1 (see the x-axis) and f2 (see the y-axis) to the frequency f3 (f1 + f2 = f3). Conversely, black indicates the gain of spectral power at frequencies f1 and f2 relative to the frequency f3. The function of f3 in the region [0,1] quantitatively describes the coupling degree between f1 and f2. When the value is closer to 1, the phase coupling degree is the largest. When the function value is 0, it indicates the smallest coupling degree.

[0059] S105. According to the spectral peak on the bispectrum diagram (judged by the heat map color), determine the frequency corresponding to this period. The reciprocal of the frequency is the millennium-scale period.

[0060] Evaluate the bispectrum of phase coupling and energy transfer between frequencies in geochemical data. The bispectrum shows that there are many phase couplings during this period, and the frequencies include but are not limited to astronomical periods. Determine the millennium-scale astronomical periods existing therein by comparing the reciprocal of the frequency corresponding to the coordinates of the reflection intensity with the time-domain index curve. For example, the reciprocals of the frequencies at the highlighted positions 0.1 and 0.3 can correspond to astronomical periods of 10 millennia and 3.33 millennia respectively.

[0061] The above-described embodiments have elaborated in detail on the technical solutions and beneficial effects of the present invention. It should be understood that the above are only specific embodiments of the present invention and are not used to limit the present invention. Any modifications, supplements, and equivalent replacements made within the scope of the principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for identifying millennium-scale cycles based on bispectral analysis, characterized in that, Including: (1) Select paleoclimate proxy indicators for detecting millennial-scale periodic climate changes, perform power spectrum analysis after preprocessing, and identify thickness cycles corresponding to astronomical orbital periods; (2) Perform Gaussian band-pass filtering on the frequencies corresponding to the identified thickness cycles, select the filtered curve to establish a time-depth conversion model, and obtain the time-domain sequence of the proxy indicator based on the time-depth conversion model; (3) Perform power spectrum analysis on the time-domain sequence to identify astronomical orbital periods and millennial-scale periods in the time-domain sequence; (4) Import the time-domain sequence into the Matlab high-order spectral analysis toolbox. Through bispectral analysis, the astronomical orbital periods in the time-domain sequence are respectively placed at frequencies f1, f2, and f3 for phase coupling. The energy is non-linearly transferred between these frequencies and redistributed on the frequency spectrum, forming a skewed and / or asymmetric cyclic geometric structure; Among them, frequencies f1 and f2 correspond to the x-axis and y-axis respectively, and f3 = f1 + f2; 2. The millennium-scale period identification method based on bispectral analysis according to claim 1, wherein (5) Compare and identify the skewed and / or asymmetric geometric structures formed in the bispectral analysis, derive the bispectral diagram of the entire time-domain curve in the time domain in step (4), use the x-axis, y-axis, and projection diagonal line to observe and identify strong response regions with large correlations, and combine the millennial-scale period identified in step (3) to further verify the existence of this millennial-scale period.

3. The millennium-scale period identification method based on bispectral analysis according to claim 2, wherein In step (1), the astronomical orbital periods refer to the long eccentricity period, short eccentricity period, slope period, and precession period, and the millennial-scale period refers to a period less than the precession period and greater than 1 millennium.

4. The method for identifying millennium-scale cycles based on bispectrum analysis according to claim 1, characterized in that In step (1), based on the ratio relationship of the eccentricity, slope, and precession periods, each corresponding thickness cycle is identified; the relationship between the thickness cycle and the corresponding depth should satisfy the estimation of the sedimentation rate of the regional geological background.

5. The millennium-scale period identification method based on bispectral analysis according to claim 1, characterized in that In step (2), the frequency corresponding to the thickness cycle is the reciprocal of the thickness.

6. The millennium-scale cycle recognition method based on bispectral analysis according to claim 1, characterized in that In step (3), identifying the astronomical orbital periods in the time-domain sequence specifically includes: the identified frequencies correspond to the long eccentricity, short eccentricity, slope, and precession periods, which are respectively: 0.00245, 0.001, 0.0025, and 0.005 cycles / millennium. In step (4), assume that the autocorrelation function of the time-domain sequence {x(n)} is expressed as: r(τ) = E{x(n)x(n + τ)} The third-order matrix formula of {x(n)} is as follows: R(τ1, τ2) = E{x(n)x(n + τ1)x(n + τ2)} B(f1,f2) = E{X(f1)X(f2)X * (f1 + f2)} where τ1 and τ2 are random times, E{} is the mathematical expectation, X(f) is the Fourier transform of {x(n)}, and X * (f1 + f2) is the complex conjugate; B(f1, f2) reflects the non-linear interaction of the astronomical signal components in the data at f1 and f2, and is considered to be the decomposition of the skewness of the signal in the frequency domain range; The bispectrum is expressed in the form of a frequency characteristic function, that is, the two-dimensional Fourier transform of R(τ1, τ2), and the formula is as follows: The bispectral coherence is used to describe the asymmetric phase coupling of the signal. The above bispectrum is normalized to obtain the bispectral coherence, as shown in the formula:

7. The method for identifying millennium-scale cycles based on bispectrum analysis according to claim 1, characterized in that In the formula, P(ω1), P(ω2), and P(ω1 + ω2) are the values of the power spectrum of x(n) at ω1, ω2, and ω1 + ω2. In step (5), if there is no millennial-scale period signal at this frequency, then the bispectral coherence coefficient is close to 0. When a millennial-scale period signal is detected, the bispectral coherence coefficient no longer remains zero. The more significant the millennial-scale period signal, the closer the bispectral coherence coefficient is to 1; Determine the frequency corresponding to this period based on the spectral peak on the bispectrum diagram, and the reciprocal of the frequency is the millennium-scale period.

8. A millennium-scale periodicity recognition device based on bispectral analysis, characterized in that, It includes a memory and one or more processors. Executable code is stored in the memory. When the one or more processors execute the executable code, it is used to implement the millennium-scale period recognition method described in any one of claims 1-7.

Images