Method for measuring shape of coal underground gasification combustion zone based on optical fiber acoustic sensing

Abstract

The present application relates to the technical field of coal underground gasification monitoring, and discloses a coal underground gasification combustion cavity shape measurement method based on optical fiber acoustic sensing, comprising: acquiring a distributed optical fiber acoustic sensing space-time data body, performing band-pass filtering processing on the space-time data body to generate a filtered acoustic signal, calculating a signal arrival time based on the filtered acoustic signal to generate an observation arrival time data set, and inverting the combustion cavity shape based on the observation arrival time data set to generate a combustion cavity boundary seismic source position coordinate set; the present application solves the problem that the three-dimensional shape of the combustion cavity is difficult to monitor in real time, and realizes high-resolution, real-time and continuous shape measurement.

Description

Technical Field

[0001] This invention relates to the field of underground coal gasification monitoring technology, and more specifically, to a method for measuring the shape of the combustion zone in underground coal gasification based on fiber optic acoustic sensing. Background Technology

[0002] Underground coal gasification is a crucial technology for converting underground coal into combustible syngas in situ. Its core process occurs in the combustion chamber of a gasifier hundreds of meters underground. The morphology, size, and expansion dynamics of this combustion chamber directly determine gasification efficiency, syngas quality, and the safety of the entire process. Accurate monitoring of the three-dimensional morphology of the combustion chamber is of great significance for optimizing the gasification process, improving syngas quality, and ensuring gasification safety.

[0003] In existing technologies, monitoring of the combustion chamber primarily relies on temperature and pressure monitoring, seismic wave CT technology, and numerical simulation methods. Temperature and pressure monitoring acquires temperature and pressure data by deploying sensors at limited locations, but it can only obtain data from a limited number of points and cannot comprehensively perceive the three-dimensional morphology of the cavity. Seismic wave CT technology uses artificially generated seismic waves and receives the seismic wave signals that penetrate the combustion chamber for tomographic imaging, but this technology is costly, cannot provide continuous real-time monitoring, and has low resolution. Numerical simulation methods simulate the evolution process of the combustion chamber by establishing geological models and gasification reaction models, but their accuracy is limited by the accuracy of the geological model and cannot provide real-time feedback.

[0004] Therefore, existing technologies have drawbacks such as limited monitoring range, high cost, inability to monitor in real time, low resolution, and limited accuracy, which leads to technical problems in the real-time, accurate, and continuous monitoring of the three-dimensional morphology of the combustion zone. Summary of the Invention

[0005] This invention provides a method for measuring the shape of the combustion zone in underground coal gasification based on fiber optic acoustic sensing, solving the technical problem of difficulty in real-time, accurate, and continuous monitoring of the three-dimensional morphology of the combustion zone in related technologies.

[0006] This invention provides a method for measuring the shape of the combustion chamber in underground coal gasification based on fiber optic acoustic sensing, comprising the following steps: Acquire a distributed fiber optic acoustic sensing spatiotemporal data volume, which contains acoustic vibration amplitude information of multiple spatial locations within a continuous time period; Bandpass filtering is applied to the spatiotemporal data volume to remove noise and interference frequency bands, while retaining the effective acoustic frequency bands generated by coal combustion at the boundary of the combustion zone, thus generating a filtered acoustic signal. The signal arrival time is calculated based on the filtered acoustic signal, and an observation arrival time dataset is generated, which contains the signal arrival times corresponding to multiple detector positions. Based on the observed arrival time dataset, the shape of the combustion zone is inverted, and a set of source location coordinates at the boundary of the combustion zone is generated. This set contains the three-dimensional coordinates of multiple source locations, which are used to describe the three-dimensional morphology of the combustion zone.

[0007] Furthermore, bandpass filtering of the spatiotemporal data volume is performed using a Butterworth filter, including the following steps: Obtain the filter performance specifications, including passband edge frequency, passband maximum attenuation, stopband edge frequency, and stopband minimum attenuation, where the passband edge frequency corresponds to the frequency range of the acoustic vibration signal generated by coal combustion at the boundary of the combustion zone. The filter order is calculated based on the attenuation requirements of the passband and stopband. The ratio of minimum stopband attenuation to maximum passband attenuation is converted into a logarithmic ratio of the ratio of stopband edge frequency to passband edge frequency through logarithmic operation. The required order near the lower cutoff frequency and the upper cutoff frequency are calculated respectively. The maximum value is taken and rounded up to obtain the final order. The normalized cutoff frequency is calculated based on the passband edge frequency and the sampling frequency, and the normalized cutoff frequency is obtained by dividing the passband edge frequency by the Nyquist frequency. The filter coefficient vector is calculated based on the filter order and the normalized cutoff frequency. The absolute values ​​of the poles of the coefficient vector are verified to be less than 1 to ensure the stability of the filter. Substituting the filter coefficients into the difference equation, the current output is calculated by a weighted combination of the current input, past inputs, and past outputs, thus generating the filtered acoustic signal.

[0008] Furthermore, the envelope thresholding method is used to calculate the signal arrival time based on the filtered acoustic signal, including the following steps: The filtered acoustic signal is discretized to obtain a discrete signal; A discrete signal is transformed by a Hilbert transform to generate an analytic signal, which consists of the original discrete signal and the imaginary part of its Hilbert transform. The envelope signal is calculated based on the analytic signal. The envelope signal is determined by the magnitude of the analytic signal and reflects the trend of the original signal amplitude changing over time. The envelope signal is smoothed by removing high-frequency fluctuations through a moving average method to generate a smooth envelope signal. The dynamic threshold is calculated based on noise statistics. The noise segment of the smooth envelope signal is selected and its mean and standard deviation are calculated. The dynamic threshold is determined by the multiple of the mean and standard deviation of the noise segment. Based on the arrival time of the signal detected by the dynamic threshold, the location of the sampling point that first exceeds the dynamic threshold is found in the smooth envelope signal. This location is multiplied by the sampling interval to obtain the signal arrival time.

[0009] Furthermore, when the envelope thresholding method cannot effectively detect the signal arrival time, the energy accumulation method is used as a backup method, including the following steps: Calculate the instantaneous energy and cumulative energy of the discrete signal. The instantaneous energy is determined by the square of the sound pressure, and the cumulative energy is the sum of the instantaneous energies from the start time to the current time. Based on the relative energy threshold, the signal arrival time is detected. The relative energy threshold is set as a preset proportion of the total energy. The sampling point position that first exceeds the relative energy threshold is found in the cumulative energy function. The signal arrival time is obtained by multiplying this position by the sampling interval. In the process of signal arrival time detection, the envelope threshold method is preferred. When there are no sampling points that make the smooth envelope greater than the dynamic threshold, the energy accumulation method is used to calculate the signal arrival time.

[0010] Furthermore, the shape of the burn-in zone is inverted based on the observed arrival time dataset using the Nelder-Mead algorithm, including the following steps: Generate a residual function, which is a weighted sum of the squares of the differences between the observed arrival time and the theoretically predicted arrival time, used to measure the degree of deviation between the current source location parameters and the actual source location; An initial simplex is generated in a multidimensional space. This initial simplex consists of multiple vertices, each of which corresponds to the coordinates of a seismic source location. Sort the vertices of the simplex according to their residual function values ​​to determine the best, second-best, second-worst, and worst points; Calculate the centroids of all vertices except the worst-case vertex; The simplex is updated through reflection, expansion, contraction, and compression operations. The appropriate operation is selected based on the relationship between the residual function values ​​of the reflection point, expansion point, and contraction point and the residual function values ​​of the current simplex vertices. Determine if the simplex satisfies the convergence condition. If it does, output the best point as the coordinates of the source location at the boundary of the combustion zone. Otherwise, return to the sorting step and continue iterating.

[0011] Furthermore, the theoretically predicted arrival time is calculated based on the detector position, the initial guessed source position, and the sound wave velocity. The theoretically predicted arrival time is equal to the distance between the detector position and the initial guessed source position divided by the sound wave velocity.

[0012] Furthermore, the reflection operation moves the worst point in the opposite direction to the centroid of the simplex. The reflection point is determined by the centroid plus the difference between the centroid and the worst point multiplied by the reflection coefficient. When the residual function value of the reflection point is less than the residual function value of the best point but greater than the residual function value of the second worst point, the worst point is replaced by the reflection point.

[0013] Furthermore, the expansion operation is triggered when the reflection point is better than the best point. The expansion point is determined by the centroid plus the difference between the reflection point and the centroid multiplied by the expansion coefficient. When the residual function value of the expansion point is less than the residual function value of the reflection point, the worst point is replaced by the expansion point; otherwise, the worst point is replaced by the reflection point.

[0014] Furthermore, the contraction operation is divided into external contraction and internal contraction. When the residual function value of the reflection point is better than that of the worst point but worse than that of the second worst point, external contraction is performed. The contraction point is determined by the centroid plus the difference between the reflection point and the centroid multiplied by the contraction coefficient. When the residual function value of the reflection point is worse than that of the worst point, internal contraction is performed. The contraction point is determined by the centroid plus the difference between the worst point and the centroid multiplied by the contraction coefficient. When the residual function value of the contraction point is less than the residual function values ​​of both the reflection point and the worst point, the worst point is replaced by the contraction point. The compression operation is triggered when the shrinking operation is ineffective. It shrinks all vertices toward the best point. The compressed vertex is determined by adding the difference between the original vertex and the best point to the best point and multiplying it by the compression coefficient.

[0015] This invention provides a coal underground gasification goaf shape measurement system based on fiber optic acoustic sensing, used to perform the aforementioned method, including: The distributed fiber optic acoustic sensor data acquisition module is used to lay fiber optic cables in underground coal gasification sites and connect them to the distributed fiber optic acoustic sensor demodulator on the ground to continuously acquire the sound wave vibration signals generated by coal combustion and generate spatiotemporal data volumes. The bandpass filter processing module is used to filter the spatiotemporal data volume, retain the effective sound wave frequency band, and generate the filtered sound wave signal. The signal arrival time calculation module is used to detect the arrival time of the signal at each detector based on the filtered acoustic signal and generate an observation arrival time dataset. The combustion zone shape inversion module is used to invert the location of the combustion zone boundary source based on the observed arrival time dataset using a derivative-free optimization algorithm, and generate a set of combustion zone boundary source location coordinates.

[0016] The beneficial effects of this invention are as follows: The method provided by this invention can perform high-resolution, real-time, and continuous shape measurement of the combustion zone in underground coal gasification, which solves the technical problem that the three-dimensional morphology of the combustion zone is difficult to monitor in real time, accurately, and continuously in the prior art, and achieves the technical effects of providing accurate decision-making basis for gasification process control, optimizing gasification efficiency, improving syngas quality, and ensuring gasification safety. Attached Figure Description

[0017] Figure 1 This is a flowchart of the method for measuring the shape of the underground coal gasification combustion zone based on fiber optic acoustic sensing according to the present invention. Figure 2 This is a block flowchart of the calculation process provided in an embodiment of the present invention; Figure 3 This is a combustion zone diagram of a microseismic source fitted to a plane section with z=0±1.2m provided in an embodiment of the present invention; Figure 4 This is a combustion zone diagram of a microseismic source fitted to a plane section with y=0±1.2m provided in an embodiment of the present invention; Figure 5 This is a combustion zone diagram of a microseismic source fitted to a plane section with x=0±1.2m provided in an embodiment of the present invention. Detailed Implementation

[0018] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, some features described in the examples may be combined in other examples.

[0019] This embodiment provides a method for measuring the shape of the underground coal gasification combustion zone based on fiber optic acoustic sensing, such as... Figure 1 As shown, the method includes the following steps: Step 100: Acquire the distributed fiber optic acoustic sensor spatiotemporal data volume; Organize the acoustic vibration signals collected by the distributed fiber optic acoustic sensing system into a spatiotemporal data volume. ,in This represents the distance channel along the optical fiber. For time, this spatiotemporal data volume contains information on the amplitude of sound wave vibrations at multiple spatial locations within a continuous time period.

[0020] Step 200: Perform bandpass filtering on the spatiotemporal data volume to generate a filtered acoustic signal; Spatiotemporal data volume The signal is input into a bandpass filter to remove noise and interference frequencies, retaining the effective acoustic frequencies generated by coal combustion at the boundary of the combustion zone, and outputting the filtered acoustic signal. .

[0021] It should be noted that bandpass filtering can be performed using a Butterworth filter. The Butterworth filter processes signals through the following steps: Step 201: Obtain the filter performance specifications; Set the passband edge frequency to and The maximum passband attenuation is Set the stopband edge frequency to and The minimum stopband attenuation is Among them, the passband edge frequency corresponds to the frequency range of acoustic vibration signals generated by coal combustion at the boundary of the combustion zone, and the stopband edge frequency corresponds to the frequency range of industrial noise and high-frequency interference that need to be suppressed.

[0022] It should be noted that the passband edge frequency and It can be set to 70Hz and 180Hz respectively, with maximum passband attenuation. It can be set to 3dB; stopband edge frequency and It can be set to 56Hz and 216Hz respectively, with minimum stopband attenuation. It can be set to 40dB.

[0023] Step 202: Calculate the filter order; Based on the attenuation requirements of the passband and stopband, the filter order is calculated using the following formula. : Specifically, the required order near the lower cutoff frequency is calculated. The required order near the upper cutoff frequency Take the maximum of the two and round it up to get the final order. .

[0024] Step 203: Calculate the normalized cutoff frequency; Based on the passband edge frequency and and sampling frequency Calculate the normalized cutoff frequency and : in, , , It is the Nyquist frequency.

[0025] Step 204: Calculate the filter coefficients; According to the filter order and normalized cutoff frequency , By converting analog filter theory into digital filter implementation, the filter coefficient vector is calculated. and .

[0026] It should be noted that, in order to ensure filter stability, it is necessary to calculate the coefficient vector. Find the poles and verify that the absolute value or magnitude of all poles is less than 1. If the absolute value or magnitude of the poles is greater than or equal to 1, the filter is unstable, and the order should be reduced. The value of is taken until the filter stabilizes.

[0027] Step 205: Generate the filtered acoustic signal using the difference equation; filter coefficients and Substituting into the difference equation, the original acoustic signal is filtered: in, This is the current output. This is the current input. It is past input. This is the past output, for a bandpass filter. , , It is the forward feedback coefficient. It is the reverse feedback coefficient, output This is the filtered sound wave signal.

[0028] Step 300: Calculate the signal arrival time based on the filtered acoustic signal and generate an observed arrival time dataset; Filtered acoustic signal The analysis identifies the arrival times of the acoustic signals at each detector, generating a dataset of observed arrival times. This dataset contains signal arrival times corresponding to multiple detector locations.

[0029] It should be noted that the signal arrival time can be calculated using the envelope thresholding method. The envelope thresholding method involves the following steps: Step 301: Discretize the filtered acoustic signal; The filtered acoustic signal Discrete signal obtained after sampling ,in The sampling interval is... It is the sampling frequency. .

[0030] Step 302: Generate an analytic signal using the Hilbert transform; For discrete signals Perform a Hilbert transform to obtain the analytic signal. ,in yes The Hilbert transform.

[0031] Step 303: Calculate the envelope signal; According to the analytical signal Calculate the envelope signal .

[0032] Step 304: Smooth the envelope signal to generate a smoothed envelope signal; The envelope signal is smoothed using a moving average method. in Output smooth envelope signal .

[0033] Step 305: Calculate the dynamic threshold based on noise statistics; Select the noise segment of the smooth envelope signal: in ; Calculate the mean of the noise segment: Sum of standard deviation: Calculate the dynamic threshold based on noise statistics. .

[0034] Step 306: Detect the signal arrival time based on the dynamic threshold; In smooth envelope signal The search for the first time exceeded the dynamic threshold sampling point location Calculate the signal arrival time .

[0035] It should be noted that when the envelope thresholding method cannot effectively detect the signal arrival time, the energy accumulation method can be used as a backup. The energy accumulation method involves the following steps: Step 307: Calculate the instantaneous energy and cumulative energy of the discrete signal; Calculate instantaneous energy ,in It's sound pressure; calculate the cumulative energy function. Calculate the total energy .

[0036] Step 308: Detect signal arrival time based on relative energy threshold; Set relative energy threshold In the cumulative energy function Find the location of the sampling point that first exceeds the relative energy threshold. Calculate the signal arrival time .

[0037] It should be noted that in the signal arrival time detection process, the envelope threshold method is preferentially used when there are no sampling points that make the envelope smooth. Greater than the threshold The energy accumulation method is used to calculate the signal arrival time.

[0038] Step 400: Based on the observed arrival time dataset, invert the shape of the combustion zone and generate a set of source location coordinates at the boundary of the combustion zone; Observation arrival time dataset The detector location information, the initial guessed source location, the sound wave velocity and other parameters are input into the derivative-free optimization algorithm. Through iterative optimization, the residual between the observed arrival time and the theoretically predicted arrival time is minimized, and the set of source location coordinates at the boundary of the combustion zone is output. This set contains the three-dimensional coordinates of multiple source locations and is used to describe the three-dimensional morphology of the combustion zone.

[0039] It should be noted that the Nelder-Mead algorithm can be used for derivative-free optimization algorithms. The Nelder-Mead algorithm is a simplex optimization algorithm that, by... Generate in 3D space by The algorithm uses a simplex composed of vertices, continuously adjusting its shape and position to gradually approach a local optimum of the objective function. The inversion process involves the following steps: Step 401: Generate the residual function; Based on the observed arrival time dataset and theoretically predicted arrival time Generate residual function: in, It is the initial guessed coordinate vector of the hypocenter location at the boundary of the combustion zone. The coordinates of each initial guess point , It is the number of data points detected. It is the first The weight of each data point It is the first Each detection arrival time, It is the first The initial guessed arrival time.

[0040] It should be noted that the theoretically predicted arrival time Calculate using the following formula: in, It is the first Each detector location It is determined by shape parameters The initial guessed location of the epicenter was determined. It is the speed of sound.

[0041] Step 402: In Generate the initial simplex in 3D space; Based on the initial guess of the epicenter location, Generate in 3D space by The initial simplex formed by the vertices , where each vertex This corresponds to the coordinates of a seismic source location.

[0042] It should be noted that the initial hypothetical source location can be determined based on the previously measured boundary location of the gasification zone. Since the growth of the underground coal gasification gasification zone is a slow and continuous process, when continuously measuring the shape of the gasification zone, the initial hypothetical source point at the gasification zone boundary can be placed on the previously measured gasification zone boundary. This ensures that the distance between the initial hypothetical point and the actual hypothetical source location is small, thus improving the accuracy of the inversion results.

[0043] Step 403: Sort the vertices of the simplex according to their residual function values; Calculate each vertex of the simplex residual function value Sort the vertices according to their residual function values ​​from smallest to largest to determine the best point. (Minimum residual function value), second best Secondary defective point and worst point (The residual function value is the largest).

[0044] Step 404: Calculate the centroid of all vertices except the worst-case vertex; Calculate the difference between the worst and worst points Other The centroid of each vertex .

[0045] Step 405: Update the simplex through reflection, expansion, contraction, and compression operations; Based on the relationship between the residual function values ​​of the reflection point, expansion point, and contraction point and the residual function value of the current simplex vertex, select the appropriate operation to update the simplex and generate a new simplex.

[0046] It should be noted that the reflection operation moves the worst-case point in the opposite direction to the simplex's centroid, exploring new directions. Reflection point. Calculate using the following formula: in, It is the centroid of all vertices except the worst-case vertex. That's the worst point. It is the reflection coefficient. When the residual function value at the reflection point... The residual function value is less than the optimal point. And greater than the residual function value of the second worst point At that time, use the reflection point Replace the worst point .

[0047] It should be noted that the expansion operation involves moving further along the reflection direction when the reflection point is better than the optimal point, accelerating the approach towards the superior region. Expansion point Calculate using the following formula: in, The reflection point has been calculated. It is the expansion coefficient. When the residual function value at the reflection point... The residual function value is less than the optimal point. When the residual function value at the expansion point is reached, the expansion operation is triggered; when the residual function value at the expansion point is reached... The residual function value less than the reflection point When using extension points Replace the worst point Otherwise, use the reflection point. Replace the worst point .

[0048] It should be noted that the contraction operation involves moving the point towards the center of gravity in small steps when the reflection point is poor, avoiding blind exploration. Based on the relationship between the reflection point and the worst-case scenario, contraction is divided into external contraction and internal contraction. When the residual function value of the reflection point... When the result is better than the worst-case scenario but worse than the second-worst-case scenario, external contraction is performed, and the contraction point is... According to the formula Calculate; the residual function value at the reflection point When the condition is worse than the worst point, internal contraction occurs; the contraction point is... According to the formula Calculation, where It is the shrinkage coefficient. The residual function value at the point of shrinkage... When the residual function value is less than that of the reflection point and the worst point, use the contraction point. Replace the worst point .

[0049] It should be noted that the compression operation, when the shrinking operation is still ineffective, shrinks all vertices towards the best point, reducing the simplex volume and focusing on a local optimum. The vertices after compression... Calculate using the following formula: in, It's the best spot. It's the compression ratio, best of all. Remain unchanged, the remaining vertices move towards shrink.

[0050] Step 406: Determine whether the simplex satisfies the convergence condition; Calculate the standard deviation or the difference between the maximum and minimum values ​​of the residual function at each vertex of the simplex, and determine if it is less than a preset convergence threshold. If the convergence condition is met, stop the iteration and output the best point. The coordinates of the source location at the boundary of the combustion zone obtained by inversion are used; if the convergence condition is not met, return to step 403 and continue iterating.

[0051] Step 407: Generate the set of source location coordinates at the boundary of the combustion zone; Steps 401 to 406 are executed for multiple seismic sources at the boundary of the combustion zone to obtain the three-dimensional coordinates of the multiple seismic source locations. These coordinates are then combined to form a set of seismic source location coordinates at the boundary of the combustion zone, which describes the three-dimensional morphology of the combustion zone.

[0052] In this embodiment of the application, in order to improve the inversion accuracy, different weights can be assigned to the observation data of different detectors in step 401 based on factors such as the distance between the detector and the seismic source and the signal quality. Data from detectors that are closer to the source and have higher signal quality are given greater weight, while data from detectors that are farther away and have lower signal quality are given less weight. This increases the sensitivity of the residual function to high-quality data and improves the accuracy of the inversion results.

[0053] In this embodiment of the application, to accelerate the convergence speed of the algorithm, the shape and size of the initial simplex can be optimized in step 402 based on prior information such as geological conditions and historical data. The vertices of the initial simplex are distributed within a reasonable range near the initially guessed epicenter location to avoid the initial simplex being too large, which would increase the number of iterations, or too small, which would lead to getting trapped in a local optimum.

[0054] In this embodiment of the application, in order to improve the stability of the algorithm, the reflection coefficient can be dynamically adjusted according to the current iteration in step 405. Expansion coefficient Shrinkage coefficient and compressibility The value of is determined by the number of coefficients. Larger coefficients are used in the early stages of iteration to accelerate the exploration process, while smaller coefficients are used in the later stages to improve convergence accuracy.

[0055] It is understood that data preprocessing methods known to those skilled in the art include data cleaning, data transformation, and data reduction. Data transformation includes type conversion and normalization and standardization. Although the dimensions and types of data were omitted in the description of the preceding embodiments, data preprocessing is a technical knowledge known to those skilled in the art and a prerequisite step in data processing. Therefore, the previously described well-known data preprocessing steps were not described independently.

[0056] The following is an example of an application of the present invention, such as Figure 2-5 As shown, the implementation process is as follows: In the underground coal gasification process, the growth of the goaf is a slow and continuous process. When using the Nelder-Mead algorithm to invert the shape of the goaf, the selection of the initial guess point is crucial for the accuracy of the calculation results. The initial guess point needs to be as close as possible to the actual hypocenter location to ensure accurate hypocenter location inversion. Since the shape of the goaf is measured in real time using fiber optic acoustic sensing, and the shape change of the goaf is slow and stable, the initial guess point for the hypocenter boundary of the goaf can be placed on the goaf boundary measured in the previous time step when continuously measuring the shape of the goaf. Because the measurement is continuous over time, the shape change of the goaf is very small, so the actual hypocenter location is very close to the initial guess point, thus ensuring the accuracy of the inversion results. Relevant calculation parameters are shown in Table 1 below.

[0057] Table 1 Calculation of relevant parameters Now assume a combustion zone is composed of a quarter sphere, half a parabola, and a bottom plane consisting of a semicircle and a parabola. Assume at time t1, this is the known shape of the combustion zone. The radius of the sphere is 10 meters, and the center coordinates of the sphere are (0, 0, 0). The equation of the sphere is: The vertex of the parabola is at (20, 0, 0), and the equation of the parabola is: Bottom plane equation: Assume that at time t2, after a period of combustion and gasification, the combustion air region expands. At this time, the radius of the combustion air sphere is 10.5 meters, the vertex coordinates of the parabola are (20.5065, 0, 0), and the bottom plane also changes accordingly, expanding by a certain point. The equation of the combustion air region at time t2 is as follows: The equation of a sphere is: The vertex of the parabola is at (20.5065, 0, 0), and the equation of the parabola is: Bottom plane equation: To ensure the detectors receive the correct signals, assume there are 200 seismic sources at the boundary of the combustion zone at time t2 (the seismic source signals are the sound waves emitted by the gasification and combustion of coal at the boundary of the combustion zone). Assume there are 40 equally spaced detection points arranged in the injection well using fiber optic acoustic sensors. Using MATLAB programming, we simulate the noisy Ricker wavelet signals emitted by the 200 seismic sources in the combustion zone at time t2 being captured by the 40 equally spaced detection points. Then each detector receives a series of signals from different seismic sources. The number of these signals is the sampling rate multiplied by the signal duration, which is 1000 in this case. Suppose we don't know the shape of the combustion zone at time t2, but we have a signal emitted by the combustion zone at time t2 that has been detected by the detection points. We can use the detected signal to invert the shape of the combustion zone at time t2. We also know the shape of the combustion zone at time t1. The shape of the combustion zone at time t1 is not significantly different from that at time t2. We set the initial guess points of 200 seismic sources on the 200 seismic sources in the combustion zone at time t1. Since the shape of the combustion zone does not change significantly, we can determine that the initial guess points are very close to the actual seismic sources, thus ensuring the accuracy of the calculation.

[0058] After obtaining the sound pressure signal detected by the detector, the signal is filtered by a Butterworth filter. Then, the signal arrival time is calculated using the envelope thresholding method (main method) and the energy accumulation method (backup method). Then, the location of the source at the boundary of the combustion zone is calculated using the Nelder-Mead algorithm, and the shape of the combustion zone is obtained through multiple locations.

[0059] The block flowchart of the entire calculation process is shown below. Figure 2 .

[0060] Since 200 uniform seismic sources were set throughout the combustion zone, the number of source points on the strictly defined x=0, y=0, z=0 cross section was relatively small. Therefore, the source points in the figure below were designed to be located at cross section x=0±1.2m, y=0±1.2m, z=0±1.2m. This increased the number of points within this range, allowing for a more complete shape of the combustion zone to be fitted. Forty detectors were used, evenly distributed along the straight line between the points (-8, 0, -4) and (8, 0, -4). The results are as follows. Figure 3 , Figure 4 and Figure 5 .

[0061] Since the microseismic signal generated at time t2 and the initial guessed points both contain some noise, it is normal for the inversion results to have some error compared to the actual shape of the combustion zone. Additionally, it should be specifically mentioned that... Figure 3Since the intercepted plane includes z=0, and z=0 is also the surface of the combustion zone, it, along with the spherical and parabolic surfaces, forms the shape of the combustion zone. Because there are 200 uniform seismic sources set throughout the combustion zone, a considerable portion of these sources are distributed to the bottom surface of the combustion zone, i.e., the z=0 plane. Therefore, this is why... Figure 3 There are many seismic sources in the plane.

[0062] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A method for measuring the shape of the combustion chamber in underground coal gasification based on fiber optic acoustic sensing, characterized in that, Includes the following steps: Acquire a distributed fiber optic acoustic sensing spatiotemporal data volume, which contains acoustic vibration amplitude information of multiple spatial locations within a continuous time period; Bandpass filtering is applied to the spatiotemporal data volume to remove noise and interference frequency bands, while retaining the effective acoustic frequency bands generated by coal combustion at the boundary of the combustion zone, thus generating a filtered acoustic signal. The signal arrival time is calculated based on the filtered acoustic signal, and an observation arrival time dataset is generated, which contains the signal arrival times corresponding to multiple detector positions. Based on the observed arrival time dataset, the shape of the combustion zone is inverted, and a set of source location coordinates at the boundary of the combustion zone is generated. This set contains the three-dimensional coordinates of multiple source locations, which are used to describe the three-dimensional morphology of the combustion zone.

2. The method according to claim 1, characterized in that, Bandpass filtering of the spatiotemporal data volume using a Butterworth filter includes the following steps: Obtain the filter performance specifications, including passband edge frequency, passband maximum attenuation, stopband edge frequency, and stopband minimum attenuation, where the passband edge frequency corresponds to the frequency range of the acoustic vibration signal generated by coal combustion at the boundary of the combustion zone. The filter order is calculated based on the attenuation requirements of the passband and stopband. The ratio of minimum stopband attenuation to maximum passband attenuation is converted into a logarithmic ratio of the ratio of stopband edge frequency to passband edge frequency through logarithmic operation. The required order near the lower cutoff frequency and the upper cutoff frequency are calculated respectively. The maximum value is taken and rounded up to obtain the final order. The normalized cutoff frequency is calculated based on the passband edge frequency and the sampling frequency, and the normalized cutoff frequency is obtained by dividing the passband edge frequency by the Nyquist frequency. The filter coefficient vector is calculated based on the filter order and the normalized cutoff frequency. The absolute values ​​of the poles of the coefficient vector are verified to be less than 1 to ensure the stability of the filter. Substituting the filter coefficients into the difference equation, the current output is calculated by a weighted combination of the current input, past inputs, and past outputs, thus generating the filtered acoustic signal.

3. The method according to claim 1, characterized in that, The envelope thresholding method is used to calculate the arrival time of the filtered acoustic signal, including the following steps: The filtered acoustic signal is discretized to obtain a discrete signal; A discrete signal is transformed by a Hilbert transform to generate an analytic signal, which consists of the original discrete signal and the imaginary part of its Hilbert transform. The envelope signal is calculated from the analytic signal. The envelope signal is determined by the magnitude of the analytic signal and reflects the trend of the original signal amplitude changing over time. The envelope signal is smoothed by removing high-frequency fluctuations through a moving average method to generate a smooth envelope signal. The dynamic threshold is calculated based on noise statistics. The noise segment of the smooth envelope signal is selected and its mean and standard deviation are calculated. The dynamic threshold is determined by the multiple of the mean and standard deviation of the noise segment. Based on the arrival time of the signal detected by the dynamic threshold, the location of the sampling point that first exceeds the dynamic threshold is found in the smooth envelope signal. This location is multiplied by the sampling interval to obtain the signal arrival time.

4. The method according to claim 3, characterized in that, When the envelope thresholding method cannot effectively detect the signal arrival time, the energy accumulation method is used as a backup method, including the following steps: Calculate the instantaneous energy and cumulative energy of the discrete signal. The instantaneous energy is determined by the square of the sound pressure, and the cumulative energy is the sum of the instantaneous energies from the start time to the current time. Based on the relative energy threshold, the signal arrival time is detected. The relative energy threshold is set as a preset proportion of the total energy. The sampling point position that first exceeds the relative energy threshold is found in the cumulative energy function. The signal arrival time is obtained by multiplying this position by the sampling interval. In the process of signal arrival time detection, the envelope threshold method is preferred. When there are no sampling points that make the smooth envelope greater than the dynamic threshold, the energy accumulation method is used to calculate the signal arrival time.

5. The method according to claim 1, characterized in that, The Nelder-Mead algorithm is used to invert the shape of the burn-out zone based on the observed arrival time dataset, including the following steps: Generate a residual function, which is a weighted sum of the squares of the differences between the observed arrival time and the theoretically predicted arrival time, used to measure the degree of deviation between the current source location parameters and the actual source location; An initial simplex is generated in a multidimensional space. This initial simplex consists of multiple vertices, each of which corresponds to the coordinates of a seismic source location. Sort the vertices of the simplex according to their residual function values ​​to determine the best, second-best, second-worst, and worst points; Calculate the centroids of all vertices except the worst-case vertex; The simplex is updated through reflection, expansion, contraction, and compression operations. The appropriate operation is selected based on the relationship between the residual function values ​​of the reflection point, expansion point, and contraction point and the residual function values ​​of the current simplex vertices. Determine if the simplex satisfies the convergence condition. If it does, output the best point as the coordinates of the source location at the boundary of the combustion zone. Otherwise, return to the sorting step and continue iterating.

6. The method according to claim 5, characterized in that, The theoretically predicted arrival time is calculated based on the detector position, the initial guessed source position, and the sound wave velocity. The theoretically predicted arrival time is equal to the distance between the detector position and the initial guessed source position divided by the sound wave velocity.

7. The method according to claim 5, characterized in that, The reflection operation moves the worst point in the opposite direction to the centroid of the simplex. The reflection point is determined by the centroid plus the difference between the centroid and the worst point multiplied by the reflection coefficient. When the residual function value of the reflection point is less than the residual function value of the best point but greater than the residual function value of the second worst point, the worst point is replaced by the reflection point.

8. The method according to claim 5, characterized in that, The expansion operation is triggered when the reflection point is better than the best point. The expansion point is determined by the centroid plus the difference between the reflection point and the centroid multiplied by the expansion coefficient. When the residual function value of the expansion point is less than the residual function value of the reflection point, the worst point is replaced by the expansion point; otherwise, the worst point is replaced by the reflection point.

9. The method according to claim 5, characterized in that, The contraction operation is divided into external contraction and internal contraction. External contraction is performed when the residual function value of the reflection point is better than that of the worst point but worse than that of the second worst point. The contraction point is determined by the centroid plus the difference between the reflection point and the centroid multiplied by the contraction coefficient. Internal contraction is performed when the residual function value of the reflection point is worse than that of the worst point. The contraction point is determined by the centroid plus the difference between the worst point and the centroid multiplied by the contraction coefficient. When the residual function value of the contraction point is less than the residual function values ​​of both the reflection point and the worst point, the worst point is replaced by the contraction point. The compression operation is triggered when the shrinking operation is ineffective. It shrinks all vertices toward the best point. The compressed vertex is determined by adding the difference between the original vertex and the best point to the best point and multiplying it by the compression coefficient.

10. A coal underground gasification goaf shape measurement system based on fiber optic acoustic sensing, used to perform the method according to any one of claims 1 to 9, characterized in that, include: The distributed fiber optic acoustic sensor data acquisition module is used to lay fiber optic cables in underground coal gasification sites and connect them to the distributed fiber optic acoustic sensor demodulator on the ground to continuously acquire the sound wave vibration signals generated by coal combustion and generate spatiotemporal data volumes. The bandpass filter processing module is used to filter the spatiotemporal data volume, retain the effective sound wave frequency band, and generate the filtered sound wave signal. The signal arrival time calculation module is used to detect the arrival time of the signal at each detector based on the filtered acoustic signal and generate an observation arrival time dataset. The combustion zone shape inversion module is used to invert the location of the combustion zone boundary source based on the observed arrival time dataset using a derivative-free optimization algorithm, and generate a set of combustion zone boundary source location coordinates.

Images