CN118518194A - A method for locating the source of intrusion events using distributed optical fiber monitoring - Google Patents
A method for locating the source of intrusion events using distributed optical fiber monitoring Download PDFInfo
- Publication number
- CN118518194A CN118518194A CN202410609038.2A CN202410609038A CN118518194A CN 118518194 A CN118518194 A CN 118518194A CN 202410609038 A CN202410609038 A CN 202410609038A CN 118518194 A CN118518194 A CN 118518194A
- Authority
- CN
- China
- Prior art keywords
- wave
- arrival
- intrusion event
- source
- positioning
- 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.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H9/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means
- G01H9/004—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means using fibre optic sensors
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/29—Graphical models, e.g. Bayesian networks
- G06F18/295—Markov models or related models, e.g. semi-Markov models; Markov random fields; Networks embedding Markov models
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
本发明提供了一种分布式光纤监测入侵事件震源定位方法，包括：读取分布式光纤调制解调信号，筛选同一入侵事件的声学信号；采用多种已知方法拾取各声学信号P波初至；建立随距离衰减的入侵事件震源定位目标函数；建立贝叶斯模型对所述入侵事件震源定位目标函数进行求解，求解过程中采用马尔可夫链蒙特卡洛方法拾取进行采样，迭代更新贝叶斯模型参数，得到自适应筛选后的P波初至和震源初定位结果；在自适应性筛选后的P波初至中添加噪音，再采用马尔可夫链蒙特卡洛方法结合贝叶斯模型迭代获取添加不确定度的震源定位结果。本发明提出的震源定位方法灵敏度高，能自适应筛选P波初至，降低单一拾取方法拾取精度不高的问题。
The present invention provides a method for locating the source of an intrusion event by monitoring it with a distributed optical fiber, comprising: reading a distributed optical fiber modulation and demodulation signal, screening the acoustic signal of the same intrusion event; using a plurality of known methods to pick up the P wave first arrival of each acoustic signal; establishing an intrusion event source location objective function that decays with distance; establishing a Bayesian model to solve the intrusion event source location objective function, using the Markov chain Monte Carlo method to pick up samples during the solution process, iteratively updating the Bayesian model parameters, and obtaining the P wave first arrival and source location results after adaptive screening; adding noise to the P wave first arrival after adaptive screening, and then using the Markov chain Monte Carlo method combined with the Bayesian model to iteratively obtain the source location result with added uncertainty. The source location method proposed by the present invention has high sensitivity, can adaptively screen the P wave first arrival, and reduces the problem of low picking accuracy of a single picking method.
Description
Technical Field
The invention relates to the field of intrusion event seismic sources, in particular to a distributed optical fiber monitoring intrusion event seismic source positioning method.
Background
The distributed optical fiber acoustic sensing technology plays an important role in the engineering fields of seismology, urban safety monitoring, slope engineering and the like due to the characteristics of high sensitivity, wide coverage range and the like. In these applications, intrusion event acoustic signal detection classification and source localization are critical. The key technology of the seismic source positioning comprises P-wave first arrival picking, seismic source positioning objective function construction and a solving method thereof.
Currently, many developments have been made in the P-wave first arrival picking technique, such as improved methods of long-short time window average ratio method, kurtosis method, red pool information criterion method, etc., waveform-based machine learning picking, etc. The methods achieve positive effects in improving the accuracy and efficiency of P-wave first arrival pickup. However, there is no method to guarantee an absolutely accurate pick-up of the first arrival of the P-wave in all cases, which may lead to failure of the source localization. Therefore, it is necessary to develop a source positioning method for adaptively screening P-wave first arrival pickup to improve accuracy of source positioning.
In the source location aspect of an intrusion event, the existing method often gives the same weight to all P-wave first arrivals. However, long-range propagating signals are affected by waveform attenuation, and the P-wave first arrivals tend to have larger errors, which means that these data should be given lower weight in source positioning. To solve this problem, it is necessary to develop a new source positioning objective function that can take into account the effect that the first arrival of the P-wave becomes worse as the propagation distance increases as a whole.
Meanwhile, the traditional seismic source positioning method often regards a certain P-wave first arrival picked up as a fixed value, and uncertainty caused by factors such as sampling frequency limitation, pickup error and the like is ignored. In practice, the picked-up P-wave first arrival should be regarded as a reference value around which a certain uncertainty interval is distributed. Therefore, the seismic source positioning method considering the uncertainty of the P-wave first arrival is developed, uncertainty of a seismic source positioning result is given, and the seismic source positioning method has important significance for defining the reliability of the seismic source positioning result.
In summary, in order to improve the performance of the distributed optical fiber acoustic sensing technology in engineering application, it is important to develop a source positioning method capable of adaptively screening the first arrival of P waves, reasonably considering the first arrival weight of P waves, and fully considering the uncertainty of the first arrival of P waves. The method not only can greatly improve the accuracy of the positioning of the seismic source, but also can strengthen the practical application value of the distributed optical fiber technology in the fields of seismology, urban safety monitoring, slope engineering and the like.
Disclosure of Invention
Aiming at the problems existing in the prior art, the distributed optical fiber monitoring intrusion event source positioning method is provided, is realized based on various P-wave first arrival pick-up and Bayesian theory, so as to improve the intrusion event source positioning precision and provide uncertainty of a source positioning result.
The technical scheme adopted by the invention is as follows: a distributed optical fiber monitoring intrusion event source positioning method comprises the following steps:
reading the distributed optical fiber modulation-demodulation signals, and screening acoustic signals of the same intrusion event;
picking up the P-wave first arrival of each acoustic signal by adopting a plurality of known methods;
establishing an intrusion event focus positioning objective function attenuated along with the distance;
Establishing a Bayesian model for the intrusion event source positioning objective function, sampling the P-wave first arrival and the source parameter of each acoustic signal picked up by each method by adopting a Markov chain Monte Carlo method, and iteratively updating the Bayesian model parameter to obtain the P-wave first arrival and the source first positioning result after self-adaptive screening;
noise is added in the P-wave first arrival after adaptive screening, and then a Markov chain Monte Carlo method is combined with Bayesian model iteration to determine a seismic source positioning result with uncertainty.
As a preferred solution, the reading the distributed optical fiber modem signal, screening the acoustic signals of the same intrusion event specifically includes:
And detecting the acoustic signals of the access intrusion event by adopting a long-short time window mean ratio method, and further screening by using a cross-correlation function of adjacent channels to obtain the acoustic signals of the same event.
As a preferred solution, the filtering using adjacent channel cross correlation functions to obtain the same event acoustic signal specifically includes:
calculating the maximum value of the cross-correlation function, wherein the cross-correlation function calculation formula is as follows:
Wherein, For the cross-correlation function value of signal x and signal y,AndTime series of signal x and signal y, respectively;
If the maximum value difference of the cross correlation functions of the adjacent signals is larger than a preset threshold value, the two acoustic signals do not belong to the same intrusion event; and vice versa, the same intrusion event.
As a preferred solution, the picking up the first arrival of each acoustic signal P-wave by using a plurality of known methods specifically includes:
And a P-wave first arrival is picked up by using an artificial pick-up method, an STA/LTA method, a high order statistics method, a red pool information criterion method, a cross correlation method, a fractal dimension method and an artificial neural network method.
As a preferred scheme, the method for establishing the intrusion event source positioning objective function decaying along with the distance specifically comprises the following steps:
Calculating the distance between the current seismic source and the straight line segment of each sensor and the propagation time of P waves;
calculating the P-wave first arrival coefficient of each sensor;
And calculating an intrusion event source positioning objective function decaying exponentially with the distance.
As a preferable scheme, the calculation method of the P-wave first arrival coefficient of each sensor is as follows:
Wherein, Is the P-wave first arrival coefficient of the kth sensor,For the number of the sensors, the number of the sensors is equal to the number of the sensors,For the variance of the propagation distance,Is the straight line segment distance of the current source from the kth sensor.
As a preferable scheme, the intrusion event focus positioning objective function decaying exponentially with distance is specifically:
Wherein, For the intrusion event source location objective function, x 0、y0、z0 is the coordinates of the source in the x, y, z directions, t 0 is the moment of occurrence of the source,For the first arrival of the P-wave for the kth acoustic signal determined by the ith pick-up method,The P-wave propagation time corresponding to the kth sensor.
As a preferable scheme, in the solving process, a Markov chain Monte Carlo method is adopted to sample the first arrival and the focus parameters of each acoustic signal P wave picked up by each method, and the iterative updating of the Bayesian model parameters specifically comprises:
Numbering model parameters, wherein the model parameters comprise P-wave first arrival and seismic source parameters of each acoustic signal picked up by each method;
Randomly generating an integer in a programming range, and randomly selecting a P-wave first arrival from acoustic signal pickup data corresponding to the integer to update if the integer is in the number range of acoustic signals; if the integer is larger than the number of acoustic signals, updating the corresponding seismic source parameters by adopting a preset updating step length and a Gaussian distribution random number with the average value of 0 and the standard deviation of 1.
As a preferable scheme, the method for obtaining the P-wave first arrival and source first positioning results after the self-adaptive screening specifically comprises the following steps:
Based on the updated model parameters, calculating a Bayesian model maximum likelihood function before and after updating the model parameters, and further calculating the receiving probability of the Bayesian model before and after updating the model parameters;
If the receiving probability is larger than or equal to the random number which satisfies 0-1 uniform distribution, the iteration Bayesian model parameter updating is successful, otherwise, the Bayesian model parameter before updating is maintained;
Repeating the sampling iteration process, taking the P-wave first arrival of the last acoustic signal selected in the iteration process as the P-wave first arrival after self-adaptive screening, and taking the average value of the last N seismic source parameters updated in the iteration process as a seismic source initial determination result.
As a preferable scheme, the noise is added in the first arrival of the P wave after the adaptive screening, and the method specifically comprises the following steps:
Taking the first arrival of the screened P wave as the mean value, and increasing the standard deviation on the basis of the mean value as To characterize the uncertainty of the first arrival of the P-wave, where a is the coefficient related to the sensor response,Distance from the current source to the straight line segment of the kth sensor.
Compared with the prior art, the beneficial effects of adopting the technical scheme are as follows: the invention utilizes the cross-correlation coefficient of adjacent channels to screen the same intrusion event signal, and has high sensitivity; the built seismic source positioning objective function of the intrusion event along with the distance attenuation can reduce the influence of the poor quality of the long-distance P wave first arrival; the developed multi-P-wave first-arrival MCMC sampling solution focus first-positioning method can adaptively screen the P-wave first-arrival, and reduces the problem of low pickup precision of a single pickup method; the provided seismic source Bayesian positioning method considering the uncertainty of the first arrival of the P wave can give out the uncertainty of the seismic source positioning result, and has good guiding significance for the reliability analysis of the seismic source positioning result.
Drawings
FIG. 1 is a flow chart of a distributed optical fiber monitoring intrusion event source positioning method provided by the invention.
Fig. 2 is an example of an acoustic signal monitored by the DAS.
Fig. 3 is a sequence of adjacent acoustic signal cross-correlation maxima.
FIG. 4 is a fiber optic sensor and intrusion event source location.
Fig. 5 is a graph of intrusion event theoretical propagation time and different method P-wave first arrival pickup times.
Fig. 6 is a plot of the sensor data positioning weights for the last iteration of the MCMC.
Fig. 7 is a sequence of iterative maximum likelihood functions for a plurality of picked source initial locations MCMC.
Fig. 8 is a sequence of MCMC iterative source positions for various picked source initial locations.
Fig. 9 is a view of the source initial positioning MCMC iterative adaptive preference time.
Fig. 10 is a P-wave first arrival noise data distribution diagram.
Fig. 11 is a sequence of iterative maximum likelihood functions for noisy data source localization MCMC.
FIG. 12 is a sequence of iterative source positions for noisy data source localization MCMC.
Fig. 13 (a), 13 (b), and 13 (c) are uncertainty of the intrusion event source positioning result in the xy, xz, yz directions, respectively.
Detailed Description
Embodiments of the present application are described in detail below, examples of which are illustrated in the accompanying drawings, wherein like or similar reference numerals refer to like or similar modules or modules having like or similar functions throughout. The embodiments described below by referring to the drawings are illustrative only and are not to be construed as limiting the application. On the contrary, the embodiments of the application include all alternatives, modifications and equivalents as may be included within the spirit and scope of the appended claims.
Because the traditional P-wave first-arrival picking method cannot ensure the accuracy of the P-wave first-arrival, the uncertainty of the P-wave first-arrival is not considered in the positioning of the seismic source, and the like. The embodiment of the invention provides a distributed optical fiber monitoring intrusion event focus positioning method based on various P-wave first-arrival pick-up and Bayesian theory, which integrates the technologies of same intrusion event screening, focus positioning objective function optimization, various P-wave first-arrival pick-up data self-adaptive screening, consideration of P-wave first-arrival uncertainty and the like, so as to improve the positioning accuracy of an intrusion event focus and provide uncertainty of a focus positioning result.
Referring to fig. 1, the method for positioning the seismic source of the intrusion event by using the distributed optical fiber monitoring comprises the following steps:
And step 1, screening acoustic signals belonging to the same intrusion event.
The signals which do not belong to the intrusion event can be effectively removed through screening the acoustic signals, so that the positioning accuracy is improved. In one embodiment, adjacent channel cross-correlation functions are used for screening, and in other embodiments, adjacent channel P-wave first arrival time differences can also be used for screening the acoustic signals of the same intrusion event. In this embodiment, the adjacent channel cross-correlation function is described as follows:
And detecting the acoustic signals of the intrusion event by using a long-short time window average ratio (STA/LTA) method for the acquired distributed optical fiber modulation-demodulation signals, and further screening by using a cross-correlation function of adjacent channels to obtain the acoustic signals of the same intrusion event. Wherein, the cross-correlation function calculation expression is:
Wherein, For the cross-correlation function value of signal x and signal y,AndTime series of signal x and signal y, respectively. At this time, the cross-correlation functions of signals 1 and 2, signals 3 and 4, …, etc. are calculated based on the sensor numbers in the system. When the maximum value of the cross correlation function of adjacent conduction is more than 5%, two acoustic signals can be considered to not belong to the same intrusion event; otherwise, the acoustic signals belonging to the same intrusion event are screened out by representing the acoustic signals belonging to the same intrusion event.
And 2, picking up the P-wave first arrival of each acoustic signal by adopting a plurality of known methods.
Because the traditional P-wave first arrival picking method cannot guarantee absolute accurate picking under all conditions, and thus the seismic source positioning failure can be caused, in this embodiment, it is proposed to pick up the P-wave first arrival of each acoustic signal by adopting a plurality of methods, and perform adaptive screening in the subsequent process. Preferably, the methods used in this embodiment include a manual pick-up method, an STA/LTA method, a higher order statistics (PAI-K) method, an erythro pool information criterion (AIC) method, a cross correlation method, a fractal dimension method, an artificial neural network method, and the like. The first arrival of the P-wave determined by the ith method for the kth acoustic signal is noted asK andThe number of acoustic signals and the method of picking up are numbered respectively.
And 3, establishing an intrusion event focus positioning objective function attenuated along with the distance.
Since the long-range propagating signal is subject to waveform-attenuated images, the P-wave first arrival tends to have a large error, and therefore, it is necessary to give these data a low weight. In the embodiment, an intrusion event source positioning ray travel time objective function which decays exponentially along with the propagation distance is established, and self-adaptive weighting of the first arrival of the P wave is realized. Wherein the objective function is calculated as follows:
assuming the current source location coordinates to be The kth sensor coordinates areK=1, 2, …, K is the number of sensors, and the straight line segment distance between the current seismic source and the kth sensor is:
The propagation time of the P wave is Wherein, the method comprises the steps of, wherein,The propagation velocity of the P wave is that the variance of the propagation distance isWherein, the method comprises the steps of, wherein,For the average value of the straight line segment distances from all sensors to the current seismic source, the P-wave first arrival coefficient of the kth sensor is
Wherein,Is the P-wave first arrival coefficient of the kth sensor,For the number of sensors, which corresponds to the number of acquired acoustic signals,For the variance of the propagation distance,Is the straight line segment distance of the current source from the kth sensor.
Thus, the intrusion event source positioning objective function decaying exponentially with distance is obtained:
。
Wherein, For the intrusion event source location objective function, x 0、y0、z0 is the coordinates of the source in the x, y, z directions, t 0 is the moment of occurrence of the source,For the first arrival of the P-wave for the kth acoustic signal determined by the ith pick-up method,The P-wave propagation time corresponding to the kth sensor.
It should be noted that in this embodiment, only the establishment of the source location objective function of the intrusion event that decays exponentially with distance is taken as an example, and in other embodiments, the establishment of the objective function may also be performed in the form of power law decay, logarithmic decay, and the like.
And 4, sampling P-wave first arrival pickup data of various methods by using the MCMC to solve the initial positioning of the seismic source.
Establishing a Bayesian model for solving an intrusion event source positioning objective function, wherein the maximum likelihood function calculation method comprises the following steps:
Wherein, Is thatThe matrix is formed by a matrix of the components,Is thatA matrix of k=1, 2, …, K; t is a matrix transpose symbol; Is that A covariance matrix is formed; Finger means Is a determinant of (2); Refers to the current model; Refers to the updated model; The kth element of (2) is 。
In this embodiment, markov Chain Monte Carlo (MCMC) is used for each model parameterX 0,y0,z0,t0), only one parameter is sampled at a time, and objective function solving is completed through sampling data and a Bayesian model. Specific: the process is as follows:
setting initial value of vibration source parameter (x 0,y0,z0,t0), and making model parameter [ ] X 0,y0,z0,t0) are numbered 1,2, …, K, k+1, k+2, k+3, k+4 in sequence; wherein K corresponds to K acoustic signals; k+1 to K+4 correspond to 4 source parameters;
Randomly generating an integer of [1, K+4], and randomly selecting one picked data from the data picked by the acoustic signals corresponding to the integer to update if the integer belongs to the integer of [1, K ]; if the number belongs to [ K+1, K+4], then use And updating the corresponding seismic source parameters, wherein p 0, mv and g are respectively the current value of the seismic source parameters, the updating step length and a Gaussian distribution random number with the average value of 0 and the standard deviation of 1. During the update process, the model parameters are updated one at a time.
In one embodiment, the update step in the xy direction is 0.5, the update step in the z direction is 0.001, and the update step at t 0 is 0.0005.
On the basis, the maximum likelihood function after updating the parameters is calculatedAnd calculates the receiving probability according to the maximum likelihood function before and after updatingComparing the random number u with the random number u satisfying 0-1 uniform distribution, ifThe iteration model parameters of this time are [ ]X 0,y0,z0,t0) is successfully updated, otherwise, the model parameters before updating are maintained. Repeating the sampling iterative process, setting the iterative times according to the requirement, and obtaining the iterative processThe last acquired value in the iteration sequence of (a) is the P-wave first arrival after self-adaptive screening, and the average value of the last 5000 x 0,y0,z0 values updated in the iteration process is used as a seismic source initial positioning result. In other embodiments, a grid search method may also be used to screen the first arrival of P waves.
Based on the primary positioning result of the seismic source, the embodiment also provides a Bayesian positioning method of the seismic source taking the uncertainty of the primary arrival of the P wave into consideration, namely:
And 5, seismic source Bayesian positioning considering the uncertainty of the first arrival of the P waves.
The P wave first arrival of self-adaptive screening is taken as the mean value, and the standard deviation is increased on the basisTo characterize the uncertainty of the first arrival of the P-wave, a being the coefficient related to the sensor response. And then, using the MCMC sampling method in the step 4, continuing to utilize the Bayesian model iterative computation x 0,y0,z0 to finish the seismic source positioning. It should be noted that, during each iteration, each P-wave first arrival is obtained from the above-mentioned noise-added data (i.e. the average value added is 0 and the standard deviation isData obtained after noise) is randomly sampled. And finally, taking the mean value of the last 5000 values of x 0,y0,z0 as a solution focus positioning result, drawing a density cloud chart for the last 5000 values of x 0,y0,z0, and representing uncertainty of the focus positioning result.
In this embodiment, after noise is added, a bayesian positioning method is still used to complete the positioning of the seismic source, and in other embodiments, an optimization method may be used to position each P-wave first arrival sample, and a plurality of initial positioning results are used to characterize uncertainty of the positioning result of the seismic source.
The following results are shown in fig. 2-13 to verify the distributed optical fiber monitoring intrusion event source positioning method provided by the invention.
Fig. 2 is an example of an acoustic signal monitored by the DAS. As can be seen, the DAS signal typically has a low signal-to-noise ratio, and the sensors 593-599, 600-612, 613-620 have similar waveforms, which can be considered as different intrusion events. The intrusion event 1 has obvious similar waveform segments (signals in a dotted line box), but the P wave first arrival is not obvious, and the P wave first arrival time difference is difficult to screen acoustic signals of the same intrusion event; the intrusion event 2 has obvious P-wave first arrival and similar waveform characteristics, which provides convenience for screening acoustic signals of the same event by using P-wave first arrival time difference and waveform cross correlation functions; the intrusion event 2 is mixed at the tail part of the intrusion event 3, so that the difficulty is increased for screening acoustic signals of the same intrusion event.
Fig. 3 is a sequence of adjacent acoustic signal cross-correlation maxima. The maximum value of the cross-correlation of signal pair number i corresponds to the maximum value of the cross-correlation sequences of the i and i +1 acoustic signals. The cross-correlation values of the acoustic signals 599, 600, 613 and 614 are obvious local minimum values, and the reduction ratio of the cross-correlation maximum value is larger than 5%, so that the acoustic signals 593-599, 600-612 and 613-620 are divided into 3 types of intrusion events, the manual analysis is consistent with that of FIG. 2, and the effectiveness of screening acoustic signals of the same event by the cross-correlation functions of adjacent channels is verified.
FIG. 4 is a fiber optic sensor and intrusion event source location. The coordinates of the intrusion event are ex=327955.64, ey=4407652.27, ez=1228.00, and the p-wave propagation velocity is 3200m/s, as shown in table 1. From the figure, the middle number sensor is closer to the intrusion event, i.e. the P-wave propagation time tends to increase when the sensor number extends to both sides.
TABLE 1 optical fiber sensor coordinates
Fig. 5 is a graph of intrusion event theoretical propagation time and different method P-wave first arrival pickup times. The theoretical propagation time curve in the graph is obtained by dividing the straight line segment distance between the intrusion event position and each sensor by the P wave propagation speed, and then the P wave first arrival is picked up by using a long-short time window average value ratio method (STA/LTA), a red pool information quantity criterion method (AIC) and a kurtosis method (PAI-K). For comparison, the P-wave first arrival pickup time is subtracted by one radix. The STA/LTA method has better P-wave first arrival pickup effect when the acoustic signal propagation distance is short, and the AIC method and the PAI-K method have larger errors due to the pickup stability problem when the acoustic signal propagation distance is short; the three pickup methods have larger P-wave first-arrival pickup errors as a whole when the acoustic signal propagation distance is longer.
Fig. 6 is a plot of the sensor data positioning weights for the last iteration of the MCMC. As can be seen from the figure, the sensor data closer to the intrusion event has significantly greater positioning weight; whereas for sensor data with a longer propagation distance, the positioning weight tends to be zero, which is compatible with the fact that the remote sensor in fig. 5 has poor quality of picked up data, the positioning of the seismic source should take a smaller weight.
Fig. 7 is a sequence of iterative maximum likelihood functions for a plurality of picked source initial locations MCMC. The maximum likelihood function is fast increased in the initial iteration stage and tends to be stable when the iteration is carried out for about 5000 times, so that the effectiveness of the established Bayesian model is shown. The maximum likelihood function sequence oscillates in the whole iteration process, and the probability of receiving the Bayesian model is caused by selecting a random number of 0-1 in each iteration, namely, a worse model with updated parameters can be received, so that convenience is provided for searching a global optimal source positioning result.
Fig. 8 is a sequence of MCMC iterative source positions for various picked source initial locations. As can be seen, the x and y coordinates vary greatly during the initial stages of the source localization iteration, while the z coordinates vary less throughout the iteration, due to the good distribution of the sensor network in the xy direction and poor distribution in the z direction. The seismic source positioning tends to be stable when the MCMC iterates for about 1000 times, and the initial positioning result is (327955.58, 4407652.30, 1228.49) m and the positioning error is 0.49m.
Fig. 9 is a view of the source initial positioning MCMC iterative adaptive preference time. As can be seen from the graph, for the sensor with a relatively close propagation distance, the self-adaptive optimized P-wave first arrival time has good consistency with the theoretical propagation time; for a sensor with a long propagation distance, the self-adaptive screening can remove P-wave first arrival pickup data with a particularly large error, but still can screen larger P-wave first arrival pickup error data. The method is characterized in that when the seismic source is positioned, the remote sensor data has weight which tends to zero, and the larger P-wave first-arrival pickup error data obtained by the self-adaptive screening has little influence on the positioning of the seismic source, so that the effectiveness of the self-adaptive screening of the P-wave first-arrival is shown.
Fig. 10 is a P-wave first arrival noise data distribution diagram. Taking the first arrival of the P wave obtained by screening as the mean value, adding Gaussian noise to the distance between the initial positioning result and the sensor and the response characteristic of the sensor, and obtaining the graph 10. From the figure, the added 20000 times Gaussian noise has good Gaussian distribution, and the noise added by the sensor data with larger propagation distance is larger overall.
Fig. 11 is a sequence of iterative maximum likelihood functions for noisy data source localization MCMC. Compared with the MCMC iterative maximum likelihood function sequence (figure 7) with the initial positioning result as input, the MCMC iterative maximum likelihood function sequence (figure 11) has a faster overall convergence speed, but the iterative sequence is integrally increased in the whole iterative process, which is caused by adding Gaussian noise to the P wave travel time at each iteration.
FIG. 12 is a sequence of iterative source positions for noisy data source localization MCMC. The figure shows that the MCMC iterative source position sequence with the noise data tends to be stable around the MCMC iteration for 2000 times due to the fact that a better iteration positioning initial value is given, but the MCMC iterative source position sequence is slightly poor in stability due to the fact that Gaussian noise is added to P wave walking time during each iteration due to the fact that the MCMC iterative source position sequence is relatively poor in stability of the MCMC iterative source position sequence with the initial positioning of the picking source. The source localization result is (327954.89, 4407652.80, 1228.34) m with a localization error of 0.98m, which is slightly larger than the initial localization error due to reduced Bayesian model convergence after the addition of Gaussian noise. The positioning error of the intrusion event is within 1m, and the positioning effect is good.
Fig. 13 (a), 13 (b), and 13 (c) are uncertainty of the intrusion event source positioning result in the xy, xz, yz directions, respectively. From the figure, the uncertainty of the seismic source positioning result in the xy direction is distributed circularly, and the uncertainty in the z direction is distributed elliptically, which is caused by the fact that the sensor network has good distribution in the xy direction and poor distribution in the z direction. The intrusion event positioning uncertainty distribution is concentrated, which indicates that the positioning result has good credibility.
In particular, according to embodiments of the present application, the processes described above with reference to flowcharts may be implemented as computer software programs. For example, embodiments of the present application include a computer program product comprising a computer program embodied on a computer readable medium, the computer program comprising program code for performing the method shown in the flowcharts.
It should be noted that, the computer readable medium shown in the embodiments of the present application may be a computer readable signal medium or a computer readable storage medium, or any combination of the two. The computer readable storage medium can be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or a combination of any of the foregoing. More specific examples of the computer-readable storage medium may include, but are not limited to: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a Random Access Memory (RAM), a read-Only Memory (ROM), an erasable programmable read-Only Memory (Erasable Programmable Read Only Memory, EPROM), a flash Memory, an optical fiber, a portable compact disc read-Only Memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device. In the present application, however, the computer-readable signal medium may include a data signal propagated in baseband or as part of a carrier wave, with the computer-readable program code embodied therein. Such a propagated data signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination of the foregoing. A computer readable signal medium may also be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device. Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to: wireless, wired, etc., or any suitable combination of the foregoing.
The flowcharts and block diagrams in the figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present application. Where each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams or flowchart illustration, and combinations of blocks in the block diagrams or flowchart illustration, can be implemented by special purpose hardware-based systems which perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.
The units involved in the embodiments of the present application may be implemented by software, or may be implemented by hardware, and the described units may also be provided in a processor. Wherein the names of the units do not constitute a limitation of the units themselves in some cases.
As another aspect, the present application also provides a computer program product or computer program comprising computer instructions stored in a computer readable storage medium. The processor of the computer device reads the computer instructions from the computer-readable storage medium and executes the computer instructions to cause the computer device to perform the distributed optical fiber monitoring intrusion event source localization method described in the above embodiments.
As another aspect, the present application also provides a computer-readable medium that may be contained in the electronic device described in the above embodiment; or may exist alone without being incorporated into the electronic device. The computer readable medium carries one or more programs which, when executed by one of the electronic devices, cause the electronic device to implement the distributed optical fiber monitoring intrusion event source localization method described in the above embodiments.
It should be noted that although in the above detailed description several modules or units of a device for action execution are mentioned, such a division is not mandatory. Indeed, the features and functions of two or more modules or units described above may be embodied in one module or unit in accordance with embodiments of the application. Conversely, the features and functions of one module or unit described above may be further divided into a plurality of modules or units to be embodied.
From the above description of embodiments, those skilled in the art will readily appreciate that the example embodiments described herein may be implemented in software, or may be implemented in software in combination with the necessary hardware. Thus, the technical solution according to the embodiments of the present application may be embodied in the form of a software product, which may be stored in a non-volatile storage medium (may be a CD-ROM, a U-disk, a mobile hard disk, etc.) or on a network, and includes several instructions to cause a computing device (may be a personal computer, a server, a touch terminal, or a network device, etc.) to perform the method according to the embodiments of the present application.
The specific meaning of the above terms in the present invention will be understood in detail by those skilled in the art; the accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the principles of the invention. The components of the embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations.
While embodiments of the present application have been shown and described above, it will be understood that the above embodiments are illustrative and not to be construed as limiting the application, and that variations, modifications, alternatives and variations may be made to the above embodiments by one of ordinary skill in the art within the scope of the application.
Claims (10)
1. The distributed optical fiber monitoring intrusion event source positioning method is characterized by comprising the following steps of:
reading the distributed optical fiber modulation-demodulation signals, and screening acoustic signals of the same intrusion event;
picking up the P-wave first arrival of each acoustic signal by adopting a plurality of known methods;
establishing an intrusion event focus positioning objective function attenuated along with the distance;
Establishing a Bayesian model to solve the intrusion event source positioning objective function, sampling the P-wave first arrival and source parameters of each acoustic signal picked up by each method by adopting a Markov chain Monte Carlo method in the solving process, and iteratively updating the Bayesian model parameters to obtain the P-wave first arrival and source initial positioning result after self-adaptive screening;
Noise is added in the P-wave first arrival after adaptive screening, and then a Markov chain Monte Carlo method is combined with a Bayesian model to obtain a seismic source positioning result with uncertainty.
2. The method for positioning a seismic source of a distributed optical fiber monitoring intrusion event according to claim 1, wherein the steps of reading the distributed optical fiber modem signals and screening acoustic signals of the same intrusion event comprise:
And detecting the acoustic signals of the access intrusion event by adopting a long-short time window mean ratio method, and further screening by using a cross-correlation function of adjacent channels to obtain the acoustic signals of the same event.
3. The method for positioning the seismic source of the intrusion event according to claim 2, wherein the filtering using the adjacent channel cross correlation function to obtain the acoustic signal of the same event specifically comprises:
calculating the maximum value of the cross-correlation function, wherein the cross-correlation function calculation formula is as follows:
Wherein, For the cross-correlation function value of signal x and signal y,AndTime series of signal x and signal y, respectively;
If the maximum value difference of the cross correlation functions of the adjacent signals is larger than a preset threshold value, the two acoustic signals do not belong to the same intrusion event; and vice versa, the same intrusion event.
4. The method for positioning the seismic source of the intrusion event according to claim 1 or 2, wherein the picking up the first arrival of each acoustic signal P-wave by using a plurality of known methods specifically comprises:
And a P-wave first arrival is picked up by using an artificial pick-up method, an STA/LTA method, a high order statistics method, a red pool information criterion method, a cross correlation method, a fractal dimension method and an artificial neural network method.
5. The method for positioning the intrusion event source according to claim 1 or 2, wherein the step of establishing the intrusion event source positioning objective function which decays along with the distance specifically comprises the steps of:
Calculating the distance between the current seismic source and the straight line segment of each sensor and the propagation time of P waves;
calculating the P-wave first arrival coefficient of each sensor;
And calculating an intrusion event source positioning objective function decaying exponentially with the distance.
6. The method for positioning the seismic source of the intrusion event according to claim 5, wherein the method for calculating the P-wave first arrival coefficient of each sensor is as follows:
Wherein, Is the P-wave first arrival coefficient of the kth sensor,For the number of the sensors, the number of the sensors is equal to the number of the sensors,For the variance of the propagation distance,Is the straight line segment distance of the current source from the kth sensor.
7. The method for positioning the intrusion event source according to claim 6, wherein the intrusion event source positioning objective function decaying exponentially with distance is specifically:
Wherein, For the intrusion event source location objective function, x 0、y0、z0 is the coordinates of the source in the x, y, z directions, t 0 is the moment of occurrence of the source,For the first arrival of the P-wave for the kth acoustic signal determined by the ith pick-up method,The P-wave propagation time corresponding to the kth sensor.
8. The method for positioning the seismic source of the intrusion event according to claim 1 or 2, wherein in the solving process, a markov chain monte carlo method is adopted to sample the first arrival and the seismic source parameters of each acoustic signal P wave picked up by each method, and the method comprises the steps of:
Numbering model parameters, wherein the model parameters comprise P-wave first arrival and seismic source parameters of all acoustic signals;
Randomly generating an integer in a programming range, and randomly selecting a P-wave first arrival from acoustic signal pickup data corresponding to the integer to update if the integer is in the number range of acoustic signals; if the integer is larger than the number of acoustic signals, updating the corresponding seismic source parameters by adopting a preset updating step length and a Gaussian distribution random number with the average value of 0 and the standard deviation of 1.
9. The method for positioning the seismic source of the intrusion event according to claim 8, wherein the step of obtaining the P-wave first arrival and the seismic source first positioning result after the adaptive screening specifically comprises the steps of:
based on the updated model parameters, calculating a Bayesian model maximum likelihood function before and after updating the model parameters, and further calculating the receiving probability of the Bayesian model before and after updating the model parameters;
If the receiving probability is larger than or equal to the random number which satisfies 0-1 uniform distribution, the iteration Bayesian model parameter updating is successful, otherwise, the Bayesian model parameter before updating is maintained;
Repeating the sampling iteration process, taking the P-wave first arrival of the last acoustic signal selected in the iteration process as the P-wave first arrival after self-adaptive screening, and taking the average value of the last N seismic source parameters updated in the iteration process as a seismic source initial determination result.
10. The method for positioning a seismic source of an intrusion event according to claim 1, wherein noise is added to the adaptively screened P-wave first arrival, specifically comprising:
taking the first arrival of the screened P wave as the mean value, and adding the mean value to be 0 and the standard deviation to be 0 on the basis To characterize the uncertainty of the first arrival of the P-wave, where a is the coefficient related to the sensor response,Distance from the current source to the straight line segment of the kth sensor.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410609038.2A CN118518194B (en) | 2024-05-16 | 2024-05-16 | A method for locating the source of intrusion events using distributed optical fiber monitoring |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410609038.2A CN118518194B (en) | 2024-05-16 | 2024-05-16 | A method for locating the source of intrusion events using distributed optical fiber monitoring |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN118518194A true CN118518194A (en) | 2024-08-20 |
| CN118518194B CN118518194B (en) | 2025-10-24 |
Family
ID=92275602
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410609038.2A Active CN118518194B (en) | 2024-05-16 | 2024-05-16 | A method for locating the source of intrusion events using distributed optical fiber monitoring |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN118518194B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN119493155A (en) * | 2024-10-29 | 2025-02-21 | 扎赉诺尔煤业有限责任公司 | A method for earthquake source location based on Bayesian time difference measurement uncertainty estimation |
Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5448531A (en) * | 1994-05-05 | 1995-09-05 | Western Atlas International | Method for attenuating coherent noise in marine seismic data |
| CN102435980A (en) * | 2011-09-15 | 2012-05-02 | 中南大学 | A Method for Localization of Acoustic Emission Source or Microseismic Source Based on Analytical Solution |
| KR101635791B1 (en) * | 2015-02-27 | 2016-07-04 | 전남대학교산학협력단 | Determination method for location and origin time of earthquake using arrival time of primary wave |
| CN110133715A (en) * | 2019-05-29 | 2019-08-16 | 长江大学 | A kind of microseism seismic source location method based on the first arrival time difference and addition of waveforms |
| CN110703313A (en) * | 2019-10-11 | 2020-01-17 | 重庆大学 | Acoustic emission event magnitude acquisition method and system considering sensor sensitivity and readable storage medium |
| CN111413733A (en) * | 2020-03-20 | 2020-07-14 | 重庆地质矿产研究院 | Mine micro-seismic positioning control system and method |
| CN111722280A (en) * | 2020-06-29 | 2020-09-29 | 重庆大学 | A Bayes location method, system and medium for acoustic emission events to remove systematic observation errors of P-wave first arrivals |
| CN111736208A (en) * | 2020-06-24 | 2020-10-02 | 重庆大学 | Bayes location method, system and medium for microseismic events based on variable weight combined P-wave and S-wave first-arrival data |
| US20210131278A1 (en) * | 2019-10-31 | 2021-05-06 | Halliburton Energy Services, Inc. | Locating passive seismic events in a wellbore using distributed acoustic sensing |
| CN116466390A (en) * | 2023-02-17 | 2023-07-21 | 南方海洋科学与工程广东省实验室(广州) | A method for real-time monitoring and positioning of earthquakes induced by large reservoirs |
| US20230288593A1 (en) * | 2022-03-14 | 2023-09-14 | Chevron U.S.A. Inc. | System and method for seismic depth uncertainty analysis |
| CN117950040A (en) * | 2022-10-20 | 2024-04-30 | 中国石油化工股份有限公司 | Seismic source positioning method, device, electronic equipment and medium based on distributed optical fiber |
-
2024
- 2024-05-16 CN CN202410609038.2A patent/CN118518194B/en active Active
Patent Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5448531A (en) * | 1994-05-05 | 1995-09-05 | Western Atlas International | Method for attenuating coherent noise in marine seismic data |
| CN102435980A (en) * | 2011-09-15 | 2012-05-02 | 中南大学 | A Method for Localization of Acoustic Emission Source or Microseismic Source Based on Analytical Solution |
| KR101635791B1 (en) * | 2015-02-27 | 2016-07-04 | 전남대학교산학협력단 | Determination method for location and origin time of earthquake using arrival time of primary wave |
| CN110133715A (en) * | 2019-05-29 | 2019-08-16 | 长江大学 | A kind of microseism seismic source location method based on the first arrival time difference and addition of waveforms |
| CN110703313A (en) * | 2019-10-11 | 2020-01-17 | 重庆大学 | Acoustic emission event magnitude acquisition method and system considering sensor sensitivity and readable storage medium |
| US20210131278A1 (en) * | 2019-10-31 | 2021-05-06 | Halliburton Energy Services, Inc. | Locating passive seismic events in a wellbore using distributed acoustic sensing |
| CN111413733A (en) * | 2020-03-20 | 2020-07-14 | 重庆地质矿产研究院 | Mine micro-seismic positioning control system and method |
| CN111736208A (en) * | 2020-06-24 | 2020-10-02 | 重庆大学 | Bayes location method, system and medium for microseismic events based on variable weight combined P-wave and S-wave first-arrival data |
| CN111722280A (en) * | 2020-06-29 | 2020-09-29 | 重庆大学 | A Bayes location method, system and medium for acoustic emission events to remove systematic observation errors of P-wave first arrivals |
| US20230288593A1 (en) * | 2022-03-14 | 2023-09-14 | Chevron U.S.A. Inc. | System and method for seismic depth uncertainty analysis |
| CN117950040A (en) * | 2022-10-20 | 2024-04-30 | 中国石油化工股份有限公司 | Seismic source positioning method, device, electronic equipment and medium based on distributed optical fiber |
| CN116466390A (en) * | 2023-02-17 | 2023-07-21 | 南方海洋科学与工程广东省实验室(广州) | A method for real-time monitoring and positioning of earthquakes induced by large reservoirs |
Non-Patent Citations (6)
| Title |
|---|
| LI XB, ET AL: "Locating single-point sources from arrival times containing large picking errors (LPEs): the virtual field optimization method (VFOM)", 《SCIENTIFIC REPORTS》, vol. 6, 31 January 2016 (2016-01-31), pages 1 - 12 * |
| WANG Y, ET AL: "Locating Mine Microseismic Events in a 3D Velocity Model through the Gaussian Beam Reverse-Time Migration Technique", 《SENSORS》, vol. 20, no. 9, 31 May 2020 (2020-05-31) * |
| ZHANG HL, ET AL: "Simultaneous Bayesian inversion for effective anisotropy parameters and source locations: a physical modelling study", 《GEOPHYSICAL JOURNAL INTERNATIONAL》, vol. 230, no. 1, 31 March 2022 (2022-03-31), pages 145 - 159 * |
| 王天韵: "地震动场实时预测方法研究", 《万方学位论文数据库》, 31 December 2018 (2018-12-31) * |
| 郑雨晴等: "联合改进作用距离和Dijkstra算法的复杂结构声发射定位方法及应用", 《黄金科学技术》, vol. 30, no. 3, 30 June 2022 (2022-06-30), pages 427 - 437 * |
| 郭宏扬: "矿山微震信号多尺度降噪及震源精确定位方法定位方法研究与应用", 《万方学位论文数据库》, 30 November 2023 (2023-11-30), pages 15 - 76 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN119493155A (en) * | 2024-10-29 | 2025-02-21 | 扎赉诺尔煤业有限责任公司 | A method for earthquake source location based on Bayesian time difference measurement uncertainty estimation |
| CN119493155B (en) * | 2024-10-29 | 2025-09-02 | 扎赉诺尔煤业有限责任公司 | A method for earthquake source location based on Bayesian uncertainty estimation of time difference measurement |
Also Published As
| Publication number | Publication date |
|---|---|
| CN118518194B (en) | 2025-10-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112949387B (en) | Intelligent anti-interference target detection method based on transfer learning | |
| CN110907991B (en) | Method, system and readable storage medium for source location based on data field potential value | |
| CN102243320A (en) | A Method for Picking Up Seismic First Arrivals | |
| CN112014883B (en) | A method, system, device and readable storage medium for microseismic source location based on Log-Cosh function | |
| CN114531729B (en) | Positioning method, system, storage medium and device based on channel state information | |
| CN118518194B (en) | A method for locating the source of intrusion events using distributed optical fiber monitoring | |
| CN113219526B (en) | A method for automatically picking and screening seismic wave first arrivals | |
| CN110703313A (en) | Acoustic emission event magnitude acquisition method and system considering sensor sensitivity and readable storage medium | |
| CN113970762A (en) | Method and system for positioning multistage interference source | |
| CA2412995A1 (en) | Seismic survey system | |
| CN120046375B (en) | A method and system for optimizing propulsion ratio in electromagnetic signal-level simulation in complex electromagnetic environments | |
| CN115932946A (en) | The training method of the magnitude estimation model, the magnitude estimation method and the magnitude estimation model | |
| CN110646846B (en) | Method, device and equipment for determining anisotropic parameters of VTI medium | |
| CN111175810A (en) | Microseismic signal arrival pickup method, device, equipment and storage medium | |
| CN115840253B (en) | Marine seismic data clock drift correction method and device | |
| CN113359161B (en) | High dynamic satellite communication signal acquisition method, device, medium and computing equipment | |
| CN114063148B (en) | A method and system for first arrival optimization of refracted waves based on Bayesian discrimination | |
| CN119493155B (en) | A method for earthquake source location based on Bayesian uncertainty estimation of time difference measurement | |
| CN117471550B (en) | Microseism event detection method and system based on Kalman filtering, storage medium and detection equipment | |
| CN118427533B (en) | Noise suppression method for extracting target frequency signal from CSAMT time series | |
| CN119475258B (en) | MEMS vibration sensor resonant frequency self-checking method, system and equipment | |
| CN119758451B (en) | Acoustic emission source localization method and system based on grid search and simulated annealing | |
| KR102839011B1 (en) | Method and system for deep learning network performance simulation in PIM architecture | |
| CN118245765B (en) | Method, device, equipment and medium for associating space-time trajectories | |
| CN114217350B (en) | Magnetic anomaly detection method, system and storage medium |
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 | ||
| TA01 | Transfer of patent application right | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20250117 Address after: No. 288, East Section of Fucheng Avenue, Chengdu High tech Zone, China (Sichuan) Pilot Free Trade Zone, Chengdu City, Sichuan Province, 610041 Applicant after: Jiangchuan Jinsha Hydropower Development Co.,Ltd. Country or region after: China Address before: No. 288, East Section of Fucheng Avenue, Chengdu High tech Zone, China (Sichuan) Pilot Free Trade Zone, Chengdu City, Sichuan Province, 610041 Applicant before: Jiangchuan Jinsha Hydropower Development Co.,Ltd. Country or region before: China Applicant before: Chongqing University |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant |