CN109581518A - One kind is without synchronizing current acquisition transmission and the acquisition of frequency domain measurement SIP data and processing method - Google Patents
One kind is without synchronizing current acquisition transmission and the acquisition of frequency domain measurement SIP data and processing method Download PDFInfo
- Publication number
- CN109581518A CN109581518A CN201910059976.9A CN201910059976A CN109581518A CN 109581518 A CN109581518 A CN 109581518A CN 201910059976 A CN201910059976 A CN 201910059976A CN 109581518 A CN109581518 A CN 109581518A
- Authority
- CN
- China
- Prior art keywords
- phase
- frequency
- current
- data
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 230000005540 biological transmission Effects 0.000 title claims abstract description 38
- 238000005259 measurement Methods 0.000 title claims abstract description 25
- 238000003672 processing method Methods 0.000 title claims abstract description 12
- 238000000034 method Methods 0.000 claims abstract description 59
- 238000004364 calculation method Methods 0.000 claims abstract description 54
- 238000010183 spectrum analysis Methods 0.000 claims abstract description 25
- 238000010606 normalization Methods 0.000 claims abstract description 24
- 230000001360 synchronised effect Effects 0.000 claims description 22
- 238000005070 sampling Methods 0.000 claims description 18
- 230000005284 excitation Effects 0.000 claims description 16
- 230000009466 transformation Effects 0.000 claims description 14
- 238000012545 processing Methods 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000004422 calculation algorithm Methods 0.000 description 16
- 238000012937 correction Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- 238000001514 detection method Methods 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 230000010287 polarization Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000010363 phase shift Effects 0.000 description 2
- 230000008054 signal transmission Effects 0.000 description 2
- 230000003595 spectral effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000000540 analysis of variance Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
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
-
- 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/02—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with propagation of electric current
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 Of Resistance Or Impedance (AREA)
Abstract
It is of the invention a kind of without synchronizing current acquisition transmission and the acquisition of frequency domain measurement SIP data and processing method, its main feature is that the initial phase and amplitude normalization coefficient of all frequency currents that can be emitted using absolute phase calculation method calculated in advance transmitter, generate initial phase list file on the basis of initial phase and amplitude normalization coefficient；Generate tranmitting frequency table and initial phase table；The phase value of the exciting current of its a certain frequency at any time then originates the acquisition time of launch time and the voltage to frequency data using it, obtains according to relative phase calculation method；The voltage amplitude and phase obtained further according to length from moving constraint all phase spectrum analysis, can access the complex resistivity at the frequency point.The present invention is acquired to high performance synchronizing current with to get rid of system and the dependence of real-time Transmission equipment or high-precision domain observations device, improves the flexibility of SIP method, realizes and smoothly carries out the observation of SIP method under the faint equal complex environments of GPS signal.
Description
Technical Field
The invention belongs to the field of geophysical electromagnetic detection, and particularly relates to a frequency domain induced polarization method (SIP method) data acquisition and processing method without synchronous current acquisition and transmission and frequency domain measurement.
Background
The frequency domain induced polarization method (SIP method) is an important branch of the induced polarization method, and the observation target is the complex resistivity of a measurement area. The method can describe the electrical structure of the underground medium of the test area more thoroughly, and is an important detailed investigation means. The SIP method can realize high-density observation in both a time domain and a frequency domain, and has the advantages of multi-parameter measurement and comprehensive information output compared with other geophysical exploration means.
At present, the existing SIP data acquisition modes and corresponding observation systems are mainly divided into two types:
(1) directly collecting the amplitude and phase of the voltage signal and the current signal, and then calculating the complex resistivity of the measuring area by using the frequency domain parameters;
(2) and synchronously acquiring time domain data of the voltage signal and the current signal, then performing Discrete Fourier Transform (DFT), and solving the complex resistivity of the measurement area by using a frequency domain parameter result of the transformation.
The first method needs to realize high-precision frequency domain measurement, and the observation system mainly comprises a high-power transmitter, a high-precision frequency domain measurement receiver, electrodes and other devices. The receiver is composed of key modules such as an imaginary component channel, a real component channel, a relevant detection unit and the like. The existing phase measurement method has high requirements on a transmitter and a receiver. Firstly, the waveform distortion of the excitation current needs to be as small as possible, the zero point needs to be stable, and meanwhile, the receiver needs to have high zero point stability. However, electromagnetic signals observed in the field are weak, and the phase measurement accuracy is difficult to guarantee. In addition, although the relative phase measurement is easy, when the environmental noise is large, the supply current needs to be increased to complete the high-precision amplitude measurement so as to improve the signal-to-noise ratio. In general, the frequency domain data acquisition method requires a relatively high accuracy and stability of the observation system, which increases the observation cost and is not conducive to field work.
The second observation method comprises a high-power transmitter, a receiver, a current detection device, a synchronous current transmission device, electrodes, a matched cable and other related peripheral equipment. The system has high requirements on instrument hardware, especially the synchronism and real-time performance of the system, the noise resistance of auxiliary devices and the like.
Disclosure of Invention
The invention aims to simplify the traditional SIP observation system, improve the flexibility and the application range of the SIP method and provide a new SIP data acquisition and processing algorithm. The dependence of the system on high-performance synchronous current acquisition and real-time transmission equipment or a high-precision frequency domain observation device is eliminated, and the observation cost is reduced; the flexibility of the SIP method is improved, and the SIP method observation can be smoothly carried out in complex environments such as weak GPS signals.
In order to achieve the purpose, the technical scheme of the invention is as follows:
a synchronous current acquisition and transmission and frequency domain measurement SIP data acquisition and processing method is characterized in that: the method comprises the steps of generating an initial phase table, calculating a relative phase, and calculating an absolute phase, wherein the relative phase calculation is phase movement related to current transmission time, and the absolute phase calculation is the initial phase of a current signal;
firstly, calculating initial phases and amplitude normalization coefficients of all frequency currents which can be transmitted by a transmitter in advance by using an absolute phase calculation method, and generating an initial phase table file on the basis of the initial phases and the amplitude normalization coefficients; generating a radio frequency table and an initial phase table;
secondly, the phase value of the excitation current with a certain frequency at any moment is obtained by utilizing the initial emission time and the acquisition time of the frequency voltage data according to a relative phase calculation method;
finally, automatically constraining the voltage amplitude and the phase obtained by the full-phase frequency spectrum analysis according to the length, namely obtaining the complex resistivity at the frequency point;
the length automatic constraint full-phase spectrum analysis comprises three parts of length automatic constraint full-phase matrix transformation, self convolution windowing and Fourier transformation.
The absolute phase calculation method comprises the following steps:
the full phase matrix is formed from the raw data and the frequency f of the excitation signalinAnd corresponding sampling rate fsAdding data length constraint N to an N-point self-convolution window, and obtaining the absolute phase of the signal by windowing Fourier transform of the N-point self-convolution window and a full phase matrix;
and the data length constraint is to constrain the length N of the full phase transformation to meet the requirement of whole period truncation. Since the excitation signal of the SIP method is mostly a series of square waves with different frequencies, N should be determined by the frequency f of the excitation signalinAnd corresponding sampling rate fsTo decide:
wherein d is1And d2Are all integers. d1The parameter d is required to satisfy the condition that the value in the symbol | | | | is an integer2The window size is guaranteed to be larger than a preset threshold σ, which is related to the sampling rate, as shown in table 1.
TABLE 1 Window Length threshold without sampling Rate
Order toThen d2Calculated according to the following formula:
where ceil () represents rounding up.
The voltage amplitude and the phase obtained by automatically constraining full-phase frequency spectrum analysis according to the length are calculated by the following steps:
(1) looking up a table for an initial phase table according to the frequency information to find a current initial phase and an amplitude normalization coefficient corresponding to the current frequency;
(2) according to the time information provided by the transmitting frequency table and the time for starting to collect the current voltage data, the time offset delta t between the current data time and the initial transmitting time of the frequency signal can be obtained, in addition, in the A/D collecting process, because the hardware starting interval can cause 128 data points to be lost, the time deviation caused by the lost data points needs to be considered when the delta t is calculated, and the value is equal to the current sampling rate fsIn connection with this, the Δ t of any time of the data segment is calculated as follows:
wherein f isinIs the current frequency;
nskipis the starting position of the currently calculated data segment in the voltage data block;
tini(fin) Is the current starting emission time of the current frequency;
tnow(fin) Is the current voltage signal start acquisition time;
(3) utilizing length automatic constraint full-phase spectrum analysis to calculate the phase and amplitude of the voltage, wherein the DFT length N is the same as the full-phase matrix of the current signal with the same frequency;
(4) the phase of the current signal, which is synchronized with the voltage signal being processed, is calculated as follows:
(5) calculating the complex resistivity phase and amplitude of the frequency point by using the formula (8) and the formula (9):
whereinAndphase values of the present voltage and current, respectively;
wherein,is the amplitude of the voltage that is applied,is the amplitude normalization coefficient, I (f)in) Is the current intensity;
and finally, correcting the calculation result to obtain the actually measured complex resistivity information.
The absolute phase solution is realized by utilizing length automatic constraint full-phase spectrum analysis.
The transmitting frequency table is a frequency-time table and comprises the following information: signal type, transmission frequency, current intensity, transmission duration, etc.
The initial phase table is an 'initial phase-amplitude-frequency table' obtained by utilizing length automatic constraint full-phase spectrum analysis calculation, and comprises the following information: signal frequency, initial phase, amplitude normalization factor, and length constraint value N.
Said phase positionThe amplitude calculation method comprises the steps of looking up an initial phase table according to frequency information, finding a current initial phase and an amplitude normalization coefficient corresponding to the current frequency, obtaining a time offset △ t between the current data time and the initial transmission time of the frequency signal according to time information provided by a transmission frequency table and the time when the current voltage data starts to be collected, and calculating the complex resistivity at a frequency point f by using the relationship between the voltage phase at the frequency point and the initial phase, frequency and time offset of the frequency currentinThe phase and amplitude of (b) can be obtained by using a normalization method.
Compared with the prior art, the invention has the substantive characteristics and remarkable effects:
the invention does not need a synchronous current acquisition and transmission and frequency domain measurement device, the phase of the current is calculated by combining an absolute phase calculation method and a relative phase calculation method with a frequency time table, and the amplitude is obtained by a normalization means. The phase calculation can obtain a high-precision phase calculation result without an additional correction means, and has the advantages of high precision and simple algorithm compared with the traditional DFT and windowed DFT. It is characterized in that:
1. the traditional SIP observation system is simplified, the dependence of the system on high-performance synchronous current acquisition and real-time transmission equipment or a high-precision frequency domain observation device is eliminated, and the observation cost is reduced;
2. the traditional complex resistivity calculation method is improved, and the phase calculation precision is improved by adopting the length automatic constraint full-phase spectrum analysis
3. The flexibility of the SIP method is improved, and SIP method observation can be smoothly carried out in complex environments such as weak GPS signals.
Drawings
Fig. 1 is an absolute phase calculation method.
Fig. 2 is a flow of data acquisition and processing by the SIP method of the present invention.
Fig. 3 illustrates an initial phase table generation method according to the present invention.
FIG. 4 is a complex resistivity calculation correction algorithm.
FIG. 5 is a comparison of calculated and theoretical complex resistivities obtained by different algorithms.
FIG. 6 is a comparison graph of calculated values and theoretical values of SIP complex resistivity under power frequency interference.
Fig. 7 shows the SIP field measurement result. The upper graph is the complex resistivity amplitude curve and the lower graph is the phase curve.
Fig. 8 is a SIP observation data processing software interface.
Detailed Description
The following detailed description of the invention refers to the accompanying drawings
Firstly, the real part and the imaginary part of the complex resistivity and the amplitude and the phase are equivalent through the analysis of the existing theory, and the real part and the imaginary part of the complex resistivity and the amplitude and the phase can be mutually converted. Accordingly, the complex resistivity of a certain frequency point can be obtained by using the voltage amplitude and phase of the frequency and the amplitude and phase of the current signal synchronized with the voltage data of the segment. In addition, the amplitude and phase of the voltage signal can be directly obtained by performing time-frequency transformation on the original sampling data. But because of the absence of synchronous current data, the spectral information of the current signal cannot be calculated in this way.
In this case, the present invention uses an absolute phase calculation method and a relative phase calculation method in combination with a frequency schedule to find the phase of the current, and the amplitude is obtained by a normalization means. The absolute phase calculation method is realized by utilizing length automatic constraint full-phase spectrum analysis.
The method can obtain the phase calculation result with high precision without additional correction means, and has the advantages of high precision and simple algorithm compared with the traditional DFT and windowed DFT.
The present invention utilizes an absolute phase calculation method, including a relative phase calculation, which refers to a phase shift associated with the current transit time, and an absolute phase calculation, which refers to the initial phase of the current signal. Firstly, calculating initial phase and amplitude normalization coefficients of all frequency currents which can be transmitted by a transmitter, and generating an initial phase table file on the basis of the initial phase normalization coefficients, namely generating a radio frequency table and an initial phase table; secondly, the phase value of the excitation current with a certain frequency at any moment is obtained by utilizing the initial emission time and the acquisition time of the frequency voltage data according to a relative phase calculation method; and finally, automatically constraining the voltage amplitude and the phase obtained by the full-phase frequency spectrum analysis according to the length to obtain the complex resistivity at the frequency point.
For a SIP transmitter, a certain frequency current x (f) is outputin) Is constant, at least for a short period of time. Thus, the phase of the signal at any time can be calculated by equation (1).
phase_It(fin)＝mod(Δt×fin)×2π+phase_Iini(fin) (1)
Wherein: f. ofinAnd phase _ Iini(fin) Respectively, the frequency and the initial phase of the current, and deltat is the time difference from the moment when the current starts to transmit to the moment when the data is collected.
The above analysis shows that once the amplitude and frequency of the transmit current are fixed, its frequency domain parameters are known at any time. In addition, the frequency domain information of the voltage signal can be obtained by time-frequency conversion using the original data. Thus, the complex resistivity parameter can be calculated according to the formula (2) and the formula (3).
phase_ρ(fin)＝phase_U(fin)-Δphase_I(fin)-phase_Iini(fin)
＝phase_U(fin)-mod(Δt×fin)×2π-phase_Iini(fin) (2)
Amp_ρ(fin)＝Amp_U(fin)/Amp_Iini(fin). (3)
Formula (2) shows that the complex resistivity is at frequency finThe phase position is composed of a voltage phase _ U at the frequency point and a starting phase _ I of the frequency currentiniFrequency finAnd a transmission time Δ t (transmission start time t)iniAnd voltage data acquisition time tnThe difference therebetween) are determined in combination. While the amplitude can be obtained using a normalization method.
Thus, without a synchronous current data acquisition and transmission device, the following parameters are used: current frequency finCurrent intensity (amplitude) A, initial phase _ IiniAnd a transmission time deltat (relative to the transmission start time), frequency domain information of the present excitation current signal can be determined. This is the theoretical basis for designing the real-time algorithm of the SIP data.
The key problem of the present invention is the calculation of phase, including relative phase calculation and absolute phase calculation. The former refers to the phase shift associated with the current transit time, while the latter refers to the initial phase of the current signal.
With stable and known signal frequency, the complex resistivity phase calculation error may come from the voltage phase error phase _ Uerr(fin) Current initial phase error phase _ Iinierr(fin) And a time error Δ terr. The first two are determined by the accuracy of the algorithm, and the third term is controlled by the accuracy of the time recording. Both the current start phase and the voltage phase are obtained using an absolute phase calculation method.
The most basic and most common phase calculation method is the fourier transform, and in the case of discrete signals, the Discrete Fourier Transform (DFT). Both the direct DFT and the windowed DFT are affected to some extent by the effects of spectral leakage. When there is no frequency error or the frequency error is small, the windowing makes the spectrum leakage effect more obvious. In addition, the window function can inhibit side lobe leakage and simultaneously cause main lobe blurring, and the effect is particularly obvious in multi-frequency signal processing with dense frequency components.
In order to obtain more accurate spectrum parameters, particularly phases, the method adopts length automatic constraint full-phase spectrum analysis, and has the advantages of small calculation amount and high phase accuracy. Because the invention has very high requirements on the phase calculation precision, some constraints are added to the full phase transformation.
The length automatic constraint full-phase spectrum analysis comprises three parts of full-phase matrix transformation, length automatic constraint self-convolution windowing and Fourier transformation.
As shown in fig. 1, the absolute phase calculation method of the present invention includes: the raw data form a full phase matrix and the frequency f of the excitation signalinAnd corresponding sampling rate fsAnd adding the data length constraint N to obtain an N-point self-convolution window, and obtaining the absolute phase of the signal by automatically constraining the N-point self-convolution window and the full phase matrix through windowed Fourier transform.
And the data length constraint is to constrain the self-convolution window length N of the full-phase transformation to meet the integral period truncation. Since the excitation signal of the SIP method is mostly a series of square waves with different frequencies, N should be determined by the frequency f of the excitation signalinAnd corresponding sampling rate fsTo decide:
wherein d is1And d2Are all integers. d1The parameter d is required to satisfy the condition that the value in the symbol | | | | is an integer2Ensuring that the window size is above a predetermined threshold sigma, which is related to the sampling rate, e.g.Shown in table 1.
TABLE 2 Window Length thresholds at No sampling Rate
Order toThen d2Calculated according to the following formula:
where ceil () represents rounding up.
The SIP method data acquisition and processing method of the invention is concretely realized by the following flows:
as shown in fig. 2, the SIP data acquisition and processing flow chart firstly generates a transmission frequency table and an initial phase table, and secondly establishes a "frequency-time table" (hereinafter referred to as "transmission frequency table"), which contains the following information: signal type, transmission frequency, current intensity, transmission time and the like. The transmission frequency table is then stored in the transmitter and the acquisition station. After the observation system starts to work, the transmitter starts to circularly transmit current in sequence according to the transmission frequency table. The receiver can then start collecting voltage data at any time. This is less limited than the current common SIP data collection method. To ensure that voltage signals of all frequencies can be acquired, the total acquisition time cannot be less than the total time of one cycle. Then, calculating the amplitude and the phase of the complex resistivity, wherein the calculating step comprises the following steps:
(1) looking up a table for an initial phase table according to the frequency information to find a current initial phase and an amplitude normalization coefficient corresponding to the current frequency;
(2) according to the time information provided by the transmitting frequency table and the time for starting to collect the current voltage dataIn addition, in the A/D acquisition process, because 128 data points are lost due to the hardware starting interval, the time deviation caused by the lost data points needs to be considered when calculating the delta t, and the value is equal to the current sampling rate fsIn connection with this, the Δ t of any time of the data segment is calculated as follows:
wherein f isinIs the current frequency;
nskipis the starting position of the currently calculated data segment in the voltage data block;
tini(fin) Is the current starting emission time of the current frequency;
tnow(fin) Is the current voltage signal start acquisition time;
(3) calculating the phase and amplitude of the voltage by utilizing length automatic constraint full-phase spectrum analysis, wherein the DFT length N is the same as the length of a full-phase transformation self-convolution window of the current signal with the same frequency;
(4) the phase of the current signal, which is synchronized with the voltage signal being processed, is calculated as follows:
(5) calculating the complex resistivity phase and amplitude of the frequency point by using the formula (8) and the formula (9):
whereinAndthe phase values of the present voltage and current, respectively.
Wherein,is the amplitude of the voltage that is applied,is the amplitude normalization coefficient, I (f)in) Is the current intensity.
And finally, correcting the calculation result to obtain the actually measured complex resistivity information.
As shown in fig. 3, the initial phase table generation method of the present invention. Before the SIP method observation is formally carried out, signals of all frequencies generated by an SIP transmitter are collected, then the initial phase and normalized amplitude coefficients of the signals are calculated by utilizing length automatic constraint full-phase spectrum analysis, and an initial phase-amplitude-frequency table (hereinafter referred to as an initial phase table) is generated on the basis of the initial phase-amplitude-frequency table, wherein the table comprises the following information: signal frequency, initial phase, amplitude normalization factor, and length constraint value N. The initial phase table file format is shown in table 2. Where the first column is the excitation current frequency and the second and third columns are the initial phase and amplitude normalization coefficients of the directly sampled data. The fourth and fifth columns are initial phase and amplitude normalization coefficients of the data processed with the harmonic filter. Because the field observation environment is complex, power frequency filtering is needed firstly when power frequency interference is strong, and a filter can cause certain influence on the amplitude and the phase of original sampling data, especially the phase, so that the frequency spectrum information after transmitting current filtering is added into the initial frequency table to meet the requirement of field real-time data processing.
TABLE 3 "initial phase Table" file format
And (5) an initial phase table generation process. The data processing procedure in the two dashed boxes is identical. Table 3 is an example of a transmit frequency table. FD represents a square wave of the signal waveform. Each frequency signal is continuously transmitted circularly. The transmission end time refers to a difference between the current frequency signal transmission end time and the first frequency signal transmission start time. For example, assuming that the 128Hz signal in Table 3 begins to be transmitted at time 00:00:00, the frequency signal ends at 00:00:50, while the 64Hz signal begins to be transmitted to 00:01:40, and so on.
TABLE 4 Transmit frequency Table
Because the measurement result of the frequency domain electromagnetic method contains the system self response, the calculation result needs to be corrected through calibration data, and therefore the real complex resistivity information of the measurement area can be finally obtained.
As shown in fig. 4, the complex resistivity calculation results are corrected. Because the frequency values used in each observation are not necessarily identical, and the calibrated frequency points are limited, the correction coefficient is calculated by using an interpolation method. According to the signal frequency finLooking up the calibration data table to find finThen calculating f by interpolationinAnd finally, correcting the original calculation result by using the calculated correction coefficient. The interpolation algorithm adopted by the invention is parabolic interpolation, so that table lookup is firstly carried out to find three and finClosest frequency point: f. ofk-1、 fkAnd fk+1. Then, the interpolation result is calculated according to the formula:
where r (f) represents the amplitude or phase calibration result at frequency f.
The complex resistivity correction formula is as follows:
φ(fin)＝φ(fin)-φc(fin)
ρ(fin)＝ρ(fin)/ρc(fin) (11)
wherein phi isc(fin) And ρc(fin) Respectively at frequency f for the observation systeminThe amplitude and phase of the self-response are calculated according to the formula (10).
The data acquisition and processing method flow is shown in fig. 1.
The comparison graph of the complex resistivity calculated value and the theoretical value can be obtained by different algorithms through simulation experiments.
As shown in fig. 5, the solid line represents the theoretical complex resistivity curve; the dotted lines represent the results obtained using the synchronous voltage and current data, and the discrete points represent the new method calculation results; "+," Δ ", and" □ "represent the results of the computation of the length auto-constrained full-phase spectral analysis, windowed DFT, and direct DFT, respectively. The upper graph is the complex resistivity amplitude and the lower graph is the phase.
In order to verify the reliability of the algorithm, simulated power frequency interference is added on the basis of the signal I to serve as excitation current. The power frequency interference consists of a set of sine waves with 50Hz and odd multiples thereof. The above simulation and data processing procedure is then repeated.
As shown in fig. 6, the complex resistivity calculation results are shown in fig. 6, and the algorithms represented by the different symbols and line types are consistent with fig. 5. The adopted data processing method is similar to that under the condition of no harmonic interference, and only power frequency filtering is added in the early stage.
The complex resistivity calculation accuracy of the above methods was discussed using an analysis of variance method. And the results are summarized in table 4. Simulation tests prove that under the condition of no power frequency interference, the calculation accuracy of the traditional SIP data processing method requiring synchronous current is higher than that of the novel algorithm provided by the invention, but under the condition of power frequency interference, the novel algorithm is superior to the traditional algorithm. For the new algorithm, the amplitude precision obtained by direct DFT, windowed DFT and length automatic constrained full-phase spectrum analysis is basically consistent, but the phase precision obtained by the length automatic constrained full-phase spectrum analysis is obviously higher than the former two.
TABLE 5 mean square error of different SIP data algorithms
As shown in fig. 7, which is a field test result of the SIP method of the present invention, the data acquisition time is about half an hour by using table 3 as a transmission frequency table, and three-channel voltage data is obtained altogether.
The complex resistivity amplitude and phase curve measured by the SIP method field test is smooth, three channels are basically consistent, and the phase of the high frequency band of the channel 1 is slightly lower than the measurement results of the other two channels.
The result shows that the SIP real-time data acquisition and processing method can smoothly complete SIP observation under the condition of no synchronous current acquisition and transmission equipment and no frequency domain measuring device. In addition, the conventional SIP method requires that voltage data and current data are strictly and synchronously acquired, and compared with the conventional SIP method, the SIP data acquisition method only acquiring the voltage data is more flexible and convenient.
As shown in fig. 8, the SIP observation data processing software interface shows the comparison between the observation results obtained by the SIP method and the full-phase spectrum analysis with the automatic length constraint.
In fig. 8, (a) and (b) are time domain waveforms and frequency domain waveforms of the original voltage sample data. (c) And (d) amplitude and phase calculations of complex resistivity. (e) Is a parameter setting and display section, where Array is a device type, "4" represents a symmetric quadrupole device; AB is the feed electrode distance; MN is the measured pole pitch. These parameters were manually entered before the start of the measurement according to the layout of the device. K is the device coefficient, which is automatically calculated by the program from the input parameters. f. ofinIs the current signal frequency, determined from the data acquisition time and the transmit frequency table, the sampling rate fsIs according to finAnd (4) determining. The solid line, the dotted line and the dashed line respectively represent the channel 1/2/3, "+" represents the calculation result of the length automatic constraint full-phase spectrum analysis of the invention, and "delta" represents the calculation result of the windowed fourier transform method.
Claims (6)
1. A synchronous current acquisition and transmission and frequency domain measurement SIP data acquisition and processing method is characterized in that: the method comprises the steps of generating an initial phase table, calculating relative phase, and calculating absolute phase, wherein the relative phase calculation is phase movement related to current transmission time, and the absolute phase calculation is the initial phase of a current signal; the SIP data acquisition and processing method comprises the following steps:
firstly, calculating initial phases and amplitude normalization coefficients of all frequency currents which can be transmitted by a transmitter in advance by using an absolute phase calculation method, and generating an initial phase table file on the basis of the initial phases and the amplitude normalization coefficients; generating a radio frequency table and an initial phase table;
secondly, the phase value of the excitation current with a certain frequency at any moment is obtained by utilizing the initial emission time and the acquisition time of the frequency voltage data according to a relative phase calculation method;
finally, automatically constraining the voltage amplitude and the phase obtained by the full-phase frequency spectrum analysis according to the length, namely obtaining the complex resistivity at the frequency point;
the length automatic constraint full-phase spectrum analysis comprises three parts of full-phase matrix transformation, length automatic constraint self-convolution windowing and Fourier transformation.
2. The method for collecting and processing SIP data without synchronous current collection and transmission and frequency domain measurement according to claim 1, wherein the method comprises the following steps: the absolute phase calculation method comprises the following steps:
the raw data form a full phase matrix and the frequency f of the excitation signalinAnd corresponding sampling rate fsAdding data length constraint N to obtain an N-point self-convolution window, and performing windowed Fourier transform on the N-point self-convolution window and a full phase matrix to obtain a signal absolute phase;
the data length constraint is to constrain the self-convolution window length N of the full phase transformation to meet the integral period truncation; since the excitation signal of the SIP method is mostly a series of square waves with different frequencies, N should be determined by the frequency f of the excitation signalinAnd corresponding sampling rate fsTo decide:
wherein d is1And d2Are all integers, d1It is necessary to satisfy the condition that the value in the symbol is integer, parameter d2To ensure that the window size is larger than a preset threshold σ, which is related to the sampling rate, as shown in table 1;
TABLE 1 Window Length threshold without sampling Rate
Order toThen d2Calculated according to the following formula:
where ceil () represents rounding up.
3. The method for collecting and processing SIP data without synchronous current collection and transmission and frequency domain measurement according to claim 1, wherein the method comprises the following steps: the voltage amplitude and the phase obtained by the length automatic constraint full-phase spectrum analysis method comprise the following calculation steps:
(1) looking up a table for an initial phase table according to the frequency information to find a current initial phase and an amplitude normalization coefficient corresponding to the current frequency;
(2) according to the time information provided by the transmitting frequency table and the time for starting to collect the current voltage data, the time offset delta t between the current data time and the initial transmitting time of the frequency signal can be obtained, in addition, in the A/D collecting process, because the hardware starting interval can cause 128 data points to be lost, the time deviation caused by the lost data points needs to be considered when the delta t is calculated, and the value is equal to the current sampling rate fsIn connection with this, the Δ t of any time of the data segment is calculated as follows:
wherein f isinIs the current frequency;
nskipis the starting position of the currently calculated data segment in the voltage data block;
tini(fin) Is the current starting emission time of the current frequency;
tnow(fin) Is the current voltage signal start acquisition time;
(3) calculating the phase and amplitude of the voltage by utilizing length automatic constraint full-phase spectrum analysis, wherein the DFT length N is the same as the length of a full-phase transformation self-convolution window of the current signal with the same frequency;
(4) the phase of the current signal, which is synchronized with the voltage signal being processed, is calculated as follows:
(5) calculating the complex resistivity phase and amplitude of the frequency point by using the formula (8) and the formula (9):
whereinAndphase values of the present voltage and current, respectively;
wherein,is the amplitude of the voltage that is applied,is the amplitude normalization coefficient, I (f)in) Is the current amperage;
and finally, correcting the calculation result to obtain the actually measured complex resistivity information.
4. The method for collecting and processing SIP data without synchronous current collection and transmission and frequency domain measurement according to claim 3, wherein the method comprises the following steps: the transmitting frequency table is a frequency-time table and comprises the following information: signal type, transmission frequency, current intensity, transmission duration, etc.
5. The method for collecting and processing SIP data without synchronous current collection and transmission and frequency domain measurement according to claim 3, wherein the method comprises the following steps: the absolute phase solution is realized by utilizing length automatic constraint full-phase spectrum analysis.
6. The method for collecting and processing SIP data without synchronous current collection and transmission and frequency domain measurement according to claim 3, wherein the method comprises the following steps: the initial phase table is an 'initial phase-amplitude-frequency table' obtained by utilizing length automatic constraint full-phase spectrum analysis calculation, and comprises the following information: signal frequency, initial phase, amplitude normalization factor, and length constraint value N.
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910059976.9A CN109581518A (en) | 2019-01-22 | 2019-01-22 | One kind is without synchronizing current acquisition transmission and the acquisition of frequency domain measurement SIP data and processing method |
| CN201911009751.9A CN110673223B (en) | 2019-01-22 | 2019-10-23 | SIP observation method without synchronous current acquisition and transmission |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910059976.9A CN109581518A (en) | 2019-01-22 | 2019-01-22 | One kind is without synchronizing current acquisition transmission and the acquisition of frequency domain measurement SIP data and processing method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109581518A true CN109581518A (en) | 2019-04-05 |
Family
ID=65917724
Family Applications (2)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910059976.9A Pending CN109581518A (en) | 2019-01-22 | 2019-01-22 | One kind is without synchronizing current acquisition transmission and the acquisition of frequency domain measurement SIP data and processing method |
| CN201911009751.9A Expired - Fee Related CN110673223B (en) | 2019-01-22 | 2019-10-23 | SIP observation method without synchronous current acquisition and transmission |
Family Applications After (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911009751.9A Expired - Fee Related CN110673223B (en) | 2019-01-22 | 2019-10-23 | SIP observation method without synchronous current acquisition and transmission |
Country Status (1)
| Country | Link |
|---|---|
| CN (2) | CN109581518A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111707714A (en) * | 2020-06-24 | 2020-09-25 | 中国地质大学（武汉） | Soil organic chloride pollution rapid investigation device and method based on complex resistivity |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115951836B (en) * | 2023-01-12 | 2023-09-05 | 上海奎芯集成电路设计有限公司 | DFI command and data channel universal controller and read-write method thereof |
Family Cites Families (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1248015C (en) * | 2003-11-25 | 2006-03-29 | 中国地质大学（北京） | Broad spectrum current drive measuring system and measuring mode |
| CN100356194C (en) * | 2005-02-01 | 2007-12-19 | 湖南继善高科技有限公司 | Chopping decoupler with exciting polarizing measurement |
| CN2781399Y (en) * | 2005-02-01 | 2006-05-17 | 湖南继善高科技有限公司 | Chopper decoupling device for induced polarization measuring |
| CN101986098B (en) * | 2010-09-21 | 2012-02-22 | 东南大学 | Tricolor projection-based Fourier transform three-dimensional measuring method |
| US10209386B2 (en) * | 2012-08-30 | 2019-02-19 | Exxonmobil Upstream Research Company | Processing methods for time division CSEM data |
| US10241101B2 (en) * | 2014-01-24 | 2019-03-26 | Schlumberger Technology Corporation | Method and apparatus for determining permittivity of rock matrix |
| CN105301610B (en) * | 2015-09-17 | 2017-12-22 | 西安空间无线电技术研究所 | A kind of Novel GPS L5 signal quick catching methods of anti-symbol saltus step |
| CN106842332B (en) * | 2016-12-25 | 2018-06-22 | 中南大学 | A kind of comprehensive electrical prospecting method |
| CN109149968B (en) * | 2018-08-27 | 2020-07-28 | 重庆西南集成电路设计有限责任公司 | Synchronous rectifier diode and synchronous rectification control circuit |
-
2019
- 2019-01-22 CN CN201910059976.9A patent/CN109581518A/en active Pending
- 2019-10-23 CN CN201911009751.9A patent/CN110673223B/en not_active Expired - Fee Related
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111707714A (en) * | 2020-06-24 | 2020-09-25 | 中国地质大学（武汉） | Soil organic chloride pollution rapid investigation device and method based on complex resistivity |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110673223A (en) | 2020-01-10 |
| CN110673223B (en) | 2021-06-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Ritzmann et al. | Comparison of measurement methods for 2–150-kHz conducted emissions in power networks | |
| US9917755B1 (en) | Providing fast radio-frequency delay measurements for envelope tracking | |
| CN101701984B (en) | Fundamental wave and harmonic wave detecting method based on three-coefficient Nuttall windowed interpolation FFT | |
| CN105137185A (en) | Frequency domain interpolation electric power harmonic wave analysis method based on discrete Fourier transform | |
| EP2406643B1 (en) | Time domain electromagnetic interference monitoring method and system | |
| CN110244116B (en) | DC instantaneous power metering circuit and quasi-synchronous calculation method thereof | |
| CN110673223B (en) | SIP observation method without synchronous current acquisition and transmission | |
| CN109407501B (en) | Time interval measuring method based on relevant signal processing | |
| CN109975771B (en) | Broadband digital channelization method based on signal third-order phase difference | |
| CN110967658B (en) | Analog input merging unit calibrator tracing method based on digital differential method | |
| WO2015154587A1 (en) | Measurement apparatus and method for ripple current | |
| WO2007068972B1 (en) | Measuring electrical impedance at various frequencies | |
| Zeng et al. | Deconvolution of sparse underwater acoustic multipath channel with a large time-delay spread | |
| CN104185271A (en) | Identification and positioning method for multiple passive intermodulation generation points | |
| CN101718816B (en) | Fundamental wave and harmonic wave detection method based on four-item coefficient Nuttall window interpolation FFT | |
| CN104849569B (en) | Dielectric loss measuring method | |
| CN101308175A (en) | Phase spectrum analyzer | |
| RU2435168C1 (en) | Method for harmonic analysis of periodic multifrequency signal | |
| CN109270404A (en) | A kind of voltage traveling wave accurate detecting method and device | |
| RU2695953C1 (en) | Method of estimating signal-to-noise ratio at input of receiving device for radio signal with digital amplitude modulation | |
| CN105319479B (en) | Two ends of electric transmission line fault localization system | |
| CN114488324B (en) | Wide-area electromagnetic method high-frequency information extraction method and system based on time domain signal reconstruction | |
| CN106603166B (en) | Vector measurement device and method for broadband modulation signal | |
| CN104849551B (en) | Harmonic phase angle analysis method | |
| CN107728101B (en) | Angular precision calibration method for microwave landing simulator |
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 | ||
| WD01 | Invention patent application deemed withdrawn after publication | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20190405 |