Advanced prediction method of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion

Abstract

The invention discloses a method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic joint tomography inversion. The adopted test monitoring arrangement scheme is closely combined with the advanced geological drilling conditions and is related to the diameter of tunnels (holes) and laneways. The implementation basis is the advanced geological drilling holes carried out in underground projects, and it is universal for various underground projects. The adopted heat source tracing method has good application prospects due to its convenient monitoring and pollution-free. The filtering method based on the combination of smooth low-pass filtering and adaptive Kalman filtering realizes the dynamic adjustment of the covariance matrix and the gain matrix of the state prediction error, has the ability to track the sudden change state, and effectively avoids the filter divergence. The method makes full use of the data information obtained by the traditional hydrogeological test, carries out the hydrogeological joint tomography inversion, can significantly reduce the multi-solution of the inverse problem, realizes the interpretation of the non-redundant information of the water-richness of the disaster-causing structure, and provides a theoretical basis for accurately exploring the disaster-causing structure ahead of the project in advance.

Description

Technical Field

[0001] The invention relates to the technical field of advanced prediction of disaster-causing structures in the field of underground engineering, and in particular to an advanced prediction method of disaster-causing structures based on heat source tracing and hydraulic combined tomography inversion. Background Art

[0002] As the construction focus of water conservancy and hydropower projects and railway and highway projects shifts to the western mountainous areas with complex terrain and geological conditions, underground engineering construction presents the characteristics of deep burial, complex geology, and frequent disasters. It is easy to encounter geological disasters such as water and mud inrush during construction. The unfavorable geological bodies developed in the rock mass have a very significant disaster-causing effect, including karst unfavorable geological bodies such as caves, karst pipes, and underground rivers, natural and artificial caves, fault fracture zones, strong weathering zones and other special rock and soil bodies. Among them, structures such as water-conducting faults and water-filled caves may become unstable under the combined action of head pressure and construction disturbance, and then induce water and mud inrush disasters, which become the main cause of water and mud inrush disasters. Underground engineering construction needs to accurately explore the water-rich disaster-causing structures in front of the project in advance to prevent them before they happen.

[0003] The traditional exploration method in the field of geophysics is one of the commonly used methods for disaster source detection. The conversion process between the parameters such as conductivity and magnetic permeability obtained by the method and the hydraulic characteristics of the formation is relatively complex and ambiguous. At the same time, the relevant conversion is often related to the specific formation conditions, which makes it difficult to provide a universal reference for specific projects. Combined with the latest research in the field of hydrogeology, the data information obtained from traditional hydrogeological tests such as tracer tests and pumping and injection tests for underground projects contains a lot of non-redundant information describing the water-richness of disaster-causing structures. Combined with the head drop and tracer travel time data, joint inversion is carried out, which provides a new idea for the advance detection of disaster-causing structures. In recent years, temperature has begun to attract attention as a new type of tracer that is cheap, convenient and pollution-free. The rapid development of double-insulated layered well drilling technology and distributed optical fiber temperature sensors has also created the possibility of high-precision, low-cost in-situ real-time measurement, providing technical support for the monitoring of head and temperature tracer data. In this context, it is urgent to combine heat source tracers with hydraulic information to study a joint inversion advance prediction method suitable for disaster-causing structures. Summary of the invention

[0004] The technical problem to be solved by the present invention is to provide a method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic joint tomography inversion, which is used to obtain the distribution of hydrogeological parameters and explore the location of water-rich disaster-causing structures, thereby providing guarantee for safe construction of the project.

[0005] In order to solve the above technical problems, the present invention provides a method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomography inversion, comprising the following steps:

[0006] (1) Based on the geological data, geophysical and drilling data of the specific engineering study area, the monitoring area in front of the tunnel face is delineated in combination with the length, depth and inclination of the advance drilling holes carried out in the advance geological exploration, targeting the water-rich structural fracture zone, karst development section and geophysical anomaly area in the short- and medium-distance forecast; then, combined with the location and number of advance drilling arrangements for underground engineering and the diameters of tunnels and lanes, the monitoring plan for the heat source tracer test and the layered source replenishment hydraulic test is determined, including the number of monitoring groups, the number of each group, and the monitoring spacing;

[0007] (2) Conduct heat source tracer tests and pumping and injection hydraulic tests according to the determined monitoring area and test monitoring plan, and monitor and record the tracer penetration curve and water head evolution curve;

[0008] (3) Perform curve smoothing and noise elimination filtering to extract tracer travel time and head drawdown data;

[0009] (4) Considering the thermal diffusion effect and dispersion effect in the temperature propagation process, the delay factor is determined based on the fastest time of the first-order derivative change of the temperature penetration time curve to achieve the correction of the tracer travel time. In order to avoid the uncertainty of the early travel time selection, the early travel time diagnostic analysis is further carried out for the travel time data set after correction;

[0010] (5) Based on the travel time dataset and the joint iterative reconstruction algorithm, iteratively reconstruct the medium characteristics, select the best reconstructed travel time and estimate the parameter distribution;

[0011] (6) Combined with the estimated hydrogeological parameter distribution obtained based on the travel time, the initial estimated mean, variance and related scale of the aquifer parameters are determined as prior information for further hydraulic tomography analysis.

[0012] (7) Combined with the simultaneous continuous linear estimation algorithm, hydraulic tomography inversion analysis is carried out, and the analysis is repeated repeatedly until the threshold setting is met. The distribution of hydrogeological parameters such as permeability coefficient, water storage coefficient and pressure conductivity is then calculated and determined. The location and distribution of water-rich disaster-causing structures are then determined in combination with the high permeability area.

[0013] Preferably, in step (1), the monitoring area is 1.2-1.5 times the depth of the longest borehole in front of the face, which is 30 to 100 meters. If the detection accuracy required by the project is high and the karst development degree and water-richness of the study area are high, the experimental monitoring spacing adopted is 0.5-1 times the hole diameter, and the maximum monitoring spacing is 2-3 times the hole diameter. If the advance drilling has revealed a potential water-rich disaster-causing structure, multiple boreholes within 3 times the hole diameter in the potential disaster-causing structure area should be combined, and the monitoring points should be encrypted as much as possible, and the monitoring spacing should be shortened to 0.3-0.5 times the minimum drilling spacing, and a new test should be carried out to conduct detailed detection to determine the scale and characteristics of the disaster-causing structure.

[0014] Preferably, in step (2), the heat source tracer test adopts the form of single-point delivery and multi-point monitoring. The hot water source delivery should be implemented in the monitoring point of the center borehole of the tunnel face. The temperature of the heat source used should be 40-60 degrees Celsius higher than the ambient background temperature. After the injection time is half of the total monitoring time, the temperature of the monitoring port closest to the injection port should rise by more than 20 degrees Celsius. After the injection of the heat source tracer test is stopped, the subsequent test can be carried out only after the temperature of the monitoring hole of the previous test is guaranteed to return to the initial background temperature. The temperature change of less than 5% is regarded as recovery.

[0015] Preferably, in step (2), the hydraulic test adopts single-point or multi-point pumping and injection, and multi-point continuous monitoring is carried out at the same time, and the experiment is carried out no less than three times; when the number of monitoring points is less than ten, the pumping and injection experiment is carried out no less than five times; when the number of monitoring points is more than ten, the number of tests is appropriately reduced; a short monitoring time interval is adopted in the early stage of the test, and the specific interval is shorter than 1% of the total time of the pumping and injection test, and as much early head data as possible is collected, and the test is stopped after the head reaches stability, and a quarter of the total time of the designed test is selected as the early stage of pumping and injection.

[0016] Preferably, in step (3), performing curve smoothing and denoising filtering to extract tracer travel time and head drawdown data specifically comprises the following steps:

[0017] (31) Research and design a smooth low-pass filter function with window weights that changes with variance, carry out pre-denoising, and further use the adaptive Kalman filter method based on conventional Kalman filtering and strong tracking Kalman filtering algorithms to adjust the covariance matrix and gain matrix of the state prediction error in real time according to the dynamic changes of noise, and carry out deep filtering and denoising;

[0018] (32) For smooth low-pass filtering, a Gaussian filter function is used to denoise the temperature tracer penetration curve. The Gaussian function variance is determined based on the relative size of the data variance in the window and the total variance of the entire data set. When the data variance in the window is greater than the total variance of the data set, the Gaussian function variance is selected as 1.5, otherwise it is selected as 0.5. Then, combined with the Gaussian function variance, the weight vector of each data in the window is determined, and the pre-denoising temperature curve result is obtained based on the weight and the measured temperature data;

[0019] (33) For smooth low-pass filtering, a fifth-order or sixth-order polynomial function is used to fit and denoise the water head curve, and the order of the polynomial used must be smaller than the window length. Then, the least squares fitting is performed on the given high-order polynomial, and the pre-denoised water head curve is obtained by stepwise moving average.

[0020] (34) Assuming that the system noise and the observation noise are uncorrelated, the system noise variance and the measurement noise variance are estimated and corrected by combining the filtering itself, and the time-varying attenuation factor l is introduced. k+1According to the dynamic changes of noise, the prediction covariance matrix U in the conventional Kalman filter is modified k+1,k , adjust the covariance matrix and gain matrix of the state prediction error in real time;

[0021]

[0022] Among them, Trace[] is to find the matrix trace, I k is the identity matrix, is the state transfer matrix, U k is the filter covariance matrix, V k-1 is the residual covariance matrix, Z k is the new information sequence, R k is the system noise covariance matrix, Q k is the measurement noise covariance matrix;

[0023] (35) Monte Carlo analysis was performed on the denoised temperature tracer penetration curve to extract the peak travel time data. Drawdown data were extracted from the denoised hydraulic head curve. One quarter of the total time after the start of the test was considered the early stage of the pumping test, and the time from the cessation of pumping to the complete recovery of the hydraulic head was considered the late stage of the pumping test. The rest of the time was considered the middle stage. Combining the early, middle and late stage data, no less than five drawdown data were selected from each group of hydraulic head curves, with at least two data selected from each of the early and late stages.

[0024] Preferably, in step (4), the peak time t of the temperature function derivative is determined B With peak travel time t D The quantitative relationship between the above two factors can be further obtained as follows:

[0025]

[0026]

[0027] t B =δt D (5)

[0028] In the above formula, and are the first-order and second-order time derivatives of the heat source temperature signal, t B and t D are the derivative peak time and peak time respectively, δ is the conversion factor, C is the convective heat transfer coefficient, D L is the thermal diffusivity;

[0029] For the travel time dataset after correction, Monte Carlo analysis was carried out to determine the specific travel time means and standard deviations, including t-5%, t-25%, t-50%, t-75% and t-100% as travel time datasets, and iterations were carried out in turn based on multiple sets of datasets.

[0030] Preferably, in step (5), the heat source temperature tracer migration process is regarded as a wavefront expansion form of a series of solute concentrations, and the initial slowness parameter field is assumed to be homogeneously distributed. According to the scale of the research domain and the inversion accuracy, the model domain is established and grid division is performed; based on the travel time data set and the joint iterative reconstruction algorithm, iterative medium feature reconstruction is performed in sequence, and the slowness correction increment of each unit is determined in combination with the ray path length of different units, and all units are weightedly corrected one by one to match the measured ray projection, and it is iterated until the travel time difference is less than a predetermined threshold, and then the slowness field is converted into a permeability parameter distribution, which is specifically converted into:

[0031]

[0032] Where R is the thermal resistivity, φ is the porosity, i is the local hydraulic gradient, and K is the permeability coefficient. Combined with the multiple travel time data sets obtained in step (6), a parameter sensitivity analysis is performed with the same number of iterations. By comparing the errors between the observed values ​​and the reconstructed values, the parameter distribution obtained by the best reconstructed travel time is selected.

[0033] Preferably, in step (5), when the number of iterations is the same, a parameter sensitivity analysis is performed, and the parameter distribution obtained by the best reconstruction travel time is selected by comparing the errors between the observed values ​​and the reconstructed values.

[0034] Preferably, in step (6), the statistical parameters of the hydrogeological parameters, including the mean and variance, are estimated based on the parameter distribution iteratively obtained by the tracer travel time, and then the parameters of the model domain are set to a non-homogeneous field that satisfies the Gaussian distribution, or the parameter distribution iteratively obtained by the tracer travel time is used as the initial field, and the statistical characteristics of the parameter field are estimated in combination with coring and geophysical exploration methods; at the same time, the relevant scales are determined based on the formation conditions obtained from the coring conditions.

[0035] Preferably, in step (7), based on the simultaneous continuous linear estimation algorithm, a geostatistical inversion analysis of the head data is performed, and the conditional effective head field is solved in combination with the initial parameter distribution to determine the mean and variance of the hydrogeological parameters, as well as the unconditional covariance and cross covariance of the parameters and the head response, and the related contribution coefficient information;

[0036] Subsequently, the first hydrogeological parameter estimation is determined, the parameter perturbation of the first iteration and the difference between the observed and simulated hydraulic heads of each group of pumping are obtained, and the conditional means of the parameters and their covariances are updated;

[0037] Based on the parameter distribution after the first iteration, the conditional effective head field is updated, the head disturbance and related covariance are updated, and the contribution coefficient of each parameter at each unit position under different pumping processes at different times is further calculated and determined. The previous parameter estimation is improved by repeated iterations until the updated parameter conditional covariance or head error meets the threshold setting.

[0038] The beneficial effects of the present invention are as follows: the implementation basis of the method is the advanced geological drilling carried out in the underground engineering, and there is no need to repeat the drilling operation, and the method is universal for underground engineering under different geological conditions; at the same time, the experimental monitoring arrangement scheme adopted is closely combined with the drilling length, depth, inclination, position and quantity conditions of the advanced geological drilling, and is related to the tunnel (hole) and the tunnel diameter, and the specific monitoring scheme can be dynamically adjusted according to the location and distribution of the geological structures such as the water-conducting fault and the water-filled karst cave revealed by the drilling; the heat source tracer method adopted has good application prospects because of its convenient monitoring and pollution-free, and The filtering method that combines smooth low-pass filtering with adaptive Kalman filtering can realize dynamic adjustment of the covariance matrix and gain matrix of the state prediction error, has the ability to track sudden change states, and effectively avoids filtering divergence; by carrying out joint tomographic inversion of hydrogeological data, it realizes the full integration of conventional data information obtained from hydrogeological tests such as tracer tests and water pumping and injection tests, significantly reduces the multi-solution nature of the inverse problem, and interprets the non-redundant information on the water-richness of the disaster-causing structures contained in the multi-source hydrogeological data, providing a theoretical basis for accurately exploring the disaster-causing structures existing in front of the project in advance. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is an overall schematic diagram of the underground engineering heat source tracing and hydraulic test monitoring arrangement of the present invention.

[0040] Figure 2 It is a schematic diagram of the decomposition of the underground engineering heat source tracing and hydraulic test monitoring arrangement of the present invention. DETAILED DESCRIPTION

[0041] A method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic joint tomography inversion includes the following steps:

[0042] (1) Based on the geological data, geophysical and drilling data of the specific engineering study area, for the water-rich structural fracture zone, karst development section and geophysical anomaly area in the short- and medium-distance forecast, combined with the length, depth and inclination of the advance drilling holes carried out by the advance geological exploration, the depth range of 1.2-1.5 times the longest borehole in front of the face is delineated as the monitoring area, usually in the range of 30 to 100 meters. Subsequently, combined with the location and number of advance drilling arrangements for underground engineering and the diameter of tunnels (holes) and roadways, the monitoring plan for heat source tracer test and layered source replenishment hydraulic test is determined, including the number of monitoring groups, the number of each group, and the monitoring spacing.

[0043] In underground projects such as tunnels, tunnels and lanes, monitoring points are arranged in the three-dimensional space in front of the project face in combination with advanced drilling construction design. Figure 1 and Figure 2 As one of the example solutions, the number of monitoring groups can be determined in combination with the number of advance drilling holes, and the specific number of monitoring points in each group can be determined according to the diameter size of tunnels, tunnels and lanes. If the project requires high detection accuracy and the karst development degree and water-richness of the study area are high, the test monitoring spacing that can be used is 0.5-1 times the hole diameter. Otherwise, the maximum monitoring spacing can be 2-3 times the hole diameter. In addition, if the advance drilling has revealed potential water-rich disaster-causing structures, multiple boreholes within 3 times the hole diameter near the potential disaster-causing structure area should be combined to increase the density of monitoring points as much as possible. The monitoring spacing can be considered to be 0.3-0.5 times the minimum drilling spacing, and re-testing and conducting detailed detection to determine the scale and characteristics of the disaster-causing structure;

[0044] (2) According to the determined monitoring area and test monitoring plan, conduct heat source tracer test and water pumping hydraulic test, and monitor and record the tracer penetration curve and head curve.

[0045] The heat source tracer test can be carried out in the form of single-point delivery and multi-point monitoring. The hot water source should be selected in the monitoring point of the hole drilled in the center of the face. The temperature of the heat source used should be 40-60 degrees Celsius higher than the ambient background temperature. After the injection time is half of the total monitoring time, the temperature of the monitoring port closest to the injection port should rise by more than 20 degrees Celsius. After stopping the injection, the subsequent test can be carried out only after the temperature of the monitoring hole after the previous test is restored to the initial background temperature. Generally, the temperature change is less than 5% and is considered to be recovered.

[0046] The hydraulic test adopts single-point or multi-point pumping and injection, and multi-point continuous monitoring at the same time. When the number of monitoring points is less than ten, the pumping and injection test shall be carried out for no less than five times; when the number of monitoring points is more than ten, the number of tests can be appropriately reduced. Since the water head changes rapidly in the early stage of the test, a short monitoring time interval should be adopted in the early stage of the test. The specific interval should be shorter than 1% of the total time of the pumping and injection test. As much early water head data as possible should be collected, and the test should be stopped after the water head reaches stability. It can be considered to select a quarter of the total time of the design test as the early stage of pumping and injection.

[0047] (3) Carry out curve smoothing and denoising filtering to extract data such as tracer travel time and head drop. Research and design a smooth low-pass filter function with window weighting that varies with variance, carry out pre-denoising, and then combine the adaptive Kalman filter method to carry out deep filtering and denoising;

[0048] For smooth low-pass filtering, the Gaussian filter function can be considered for denoising the temperature tracer penetration curve. In order to ensure the elimination of noise while retaining as many curve features as possible, the Gaussian function variance can be determined based on the relative size of the data variance in the window and the total variance of the overall data set. When the data variance in the window is greater than the total variance, the Gaussian function variance is selected as 1.5, otherwise it is taken as 0.5. Then, combined with the Gaussian function variance, the weight vector of each data in the window is determined, and the pre-denoising result is obtained based on the weight and the measured temperature data; on the other hand, the fifth-order or sixth-order polynomial function can be considered for fitting denoising for the head curve, and the order of the polynomial used must be less than the window length. The weighting coefficient is determined based on the proportion of the data variance in the window to the total variance of the overall data set, and then the least squares fitting is performed on the given high-order polynomial, and the pre-denoising head curve is obtained by step-by-step moving average.

[0049] For adaptive Kalman filtering, it is assumed that the system noise and the observation noise are uncorrelated. The system noise variance and the measurement noise variance are estimated and corrected by combining the filter itself. The time-varying attenuation factor l is introduced. k+1 According to the dynamic changes of noise, the prediction covariance matrix U in the conventional Kalman filter is modified k+1,k , adjust the covariance matrix and gain matrix of the state prediction error in real time to ensure that the gain matrix can be automatically adjusted according to the uncertainty of the model and the changes in external noise, so that the filter has the ability to track the sudden change state and avoid filter divergence.

[0050]

[0051]

[0052] U k is the filter covariance matrix, is the state transfer matrix, Q k is the covariance matrix of the measurement noise, I is the identity matrix, Trace[] is the matrix trace, Z k is the new information sequence, R k is the system noise covariance matrix, V k-1 is the residual covariance matrix.

[0053] Monte Carlo analysis was performed to extract peak travel time data from the temperature tracer penetration curve after noise elimination. Drawdown data were extracted from the head curve after noise elimination. Generally, one quarter of the total time after the start of the test was considered as the early stage of the pumping test, and the period from the cessation of pumping to the complete recovery of the head was considered as the late stage of the pumping test, and the rest of the time was considered as the middle stage. Combining the early, middle and late stage data, no less than five drawdown data were selected in each group of head curves, of which at least two data were selected in the early and late stages.

[0054] (4) Considering the thermal diffusion effect and dispersion effect in the temperature propagation process, the delay factor is determined based on the fastest time of the first-order derivative change of the temperature penetration time curve to achieve the correction of the tracer travel time. In order to avoid the uncertainty of the early travel time selection, the early travel time diagnostic analysis is further carried out on the travel time data set after correction.

[0055] According to the heat convection diffusion control equation, the quantitative relationship between the peak time of the temperature function derivative and the peak travel time is determined, and the conversion factor can be further obtained as follows:

[0056]

[0057]

[0058] t B =δt D (5)

[0059] In the above formula, and are the first-order and second-order time derivatives of the heat source temperature signal, t B and t D are the derivative peak time and peak time respectively, δ is the conversion factor, C is the convective heat transfer coefficient, D L is the thermal diffusivity.

[0060] For the travel time dataset after correction, Monte Carlo analysis was carried out to determine the specific travel time means and standard deviations, including t-5%, t-25%, t-50%, t-75% and t-100% as travel time datasets, and iterations were carried out in turn based on multiple sets of datasets.

[0061] (5) Based on the travel time dataset and the joint iterative reconstruction algorithm, the medium characteristics are iteratively reconstructed and the optimal reconstructed travel time and its estimated parameter distribution are selected.

[0062] The specific approach is to regard the heat source temperature tracing migration process as a series of solute concentration wavefront expansion forms, assuming that the initial slowness parameter field is homogeneously distributed, and establish the model domain and perform grid division according to the scale of the research domain and the inversion accuracy.

[0063] Subsequently, based on the temperature travel time data set and the joint iterative reconstruction algorithm, the medium characteristics are reconstructed iteratively. Combined with the ray path lengths of different units, the slowness correction increment of each unit is determined, and all units are weighted corrected one by one to match the measured ray projections. Iterations are performed until the travel time difference is less than the predetermined threshold, and then the slowness field is converted into the permeability parameter distribution. The specific conversion can be referred to:

[0064]

[0065] Among them, R is the thermal resistance coefficient, φ is the porosity, i is the local hydraulic gradient, and K is the permeability coefficient. When the number of iterations is the same, a parameter sensitivity analysis is carried out, and the parameter distribution obtained by the best reconstruction travel time is selected by comparing the errors between the observed values ​​and the reconstructed values.

[0066] (6) Based on the estimated hydrogeological parameter distribution obtained during the trip, determine the initial estimated mean, variance, and related scale information of the aquifer parameters.

[0067] In some implementation cases, the statistical parameters of the hydrogeological parameters, including the mean and variance, are estimated based on the parameter distribution obtained by iterative tracer travel time, and then the parameters of the model domain are set to a non-homogeneous field that satisfies the Gaussian distribution; or the parameter distribution obtained by iterative tracer travel time can be used as the initial field, and combined with drilling coring and geophysical exploration, the statistical characteristics of the parameter field can be estimated. In addition, the determination of the relevant scale should fully consider the formation conditions obtained from the drilling coring situation.

[0068] (7) Combined with the simultaneous continuous linear estimation algorithm, the inversion analysis of hydrogeological parameters is carried out, and the parameters are iterated repeatedly until the threshold setting is met. Then, the distribution of parameters such as permeability coefficient, water storage coefficient and pressure conductivity is calculated and determined, and the location and distribution of water-rich disaster-causing structures are delineated in combination with the high permeability area.

[0069] Specifically, based on the simultaneous continuous linear estimation algorithm, geostatistical inversion analysis of the head data is performed, and combined with the initial parameter distribution, the conditional effective head field is solved to determine the mean and variance of the hydrogeological parameters, as well as the unconditional covariance and cross-covariance of its parameters and the head response, the correlation contribution coefficient and other information.

[0070] Subsequently, the first hydrogeological parameter estimation is determined, the parameter perturbation of the first iteration and the difference between the observed head and the simulated head of each group of pumping are obtained, and the parameter conditional mean and its covariance are updated. Based on the parameter distribution after the first iteration, the conditional effective head field is updated, the head perturbation and related covariance are updated, and the contribution coefficient of the head perturbation to each parameter at each unit position under different pumping processes at different times is further calculated and determined, and the previous parameter estimation is improved by repeated iterations until the updated parameter conditional covariance or head error meets the threshold setting.

Claims

1. A method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic joint tomography inversion, characterized in that: The steps include: (1) Based on the geological data, geophysical and drilling data of the specific engineering study area, the monitoring area in front of the tunnel face is delineated in combination with the length, depth and inclination of the advance drilling hole carried out in advance geological exploration, targeting the water-rich structural fracture zone, karst development section and geophysical anomaly areas in the short- and medium-distance forecast; Subsequently, combined with the location and quantity of the advance drilling arrangement of underground projects and the diameters of tunnels and laneways, the monitoring plan for the heat source tracer test and the layered source replenishment hydraulic test was determined, including the number of monitoring groups, the number of each group, and the monitoring spacing; (2) Conduct heat source tracer tests and pumping and injection hydraulic tests according to the determined monitoring area and test monitoring plan, and monitor and record the tracer penetration curve and water head evolution curve; (3) Perform curve smoothing and noise elimination filtering to extract tracer travel time and head drawdown data; (4) Considering the thermal diffusion effect and dispersion effect in the temperature propagation process, the delay factor is determined based on the fastest time of the first-order derivative change of the temperature penetration time curve to achieve the correction of the tracer travel time; in order to avoid the uncertainty of the early travel time selection, the early travel time diagnostic analysis is further carried out on the travel time data set after correction; (5) Based on the travel time dataset and the joint iterative reconstruction algorithm, iteratively reconstruct the medium characteristics, select the best reconstructed travel time and estimate the parameter distribution; (6) Based on the estimated hydrogeological parameter distribution obtained during the trip, the initial estimated mean, variance and related scale of the aquifer parameters are determined as prior information for further hydraulic tomography analysis; (7) Combined with the simultaneous continuous linear estimation algorithm, hydraulic tomography inversion analysis is carried out, and the analysis is repeated repeatedly until the threshold setting is met. The distribution of hydrogeological parameters such as permeability coefficient, water storage coefficient and pressure conductivity is then calculated and determined. The location and distribution of water-rich disaster-causing structures are then determined in combination with the high permeability area.

2. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1 is characterized in that: In step (1), the monitoring area is 1.2-1.5 times the depth of the longest borehole in front of the face, which is 30 to 100 meters. If the project requires high detection accuracy and the study area has a high degree of karst development and strong water-richness, the test monitoring spacing used is 0.5-1 times the hole diameter, and the maximum monitoring spacing is 2-3 times the hole diameter. If the advance drilling has revealed a potential water-rich disaster-causing structure, multiple boreholes within 3 times the hole diameter in the potential disaster-causing structure area should be combined, and the monitoring points should be encrypted as much as possible, and the monitoring spacing should be shortened to 0.3-0.5 times the minimum drilling spacing. Re-test and conduct detailed detection to determine the scale and characteristics of the disaster-causing structure.

3. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1 is characterized in that: In step (2), the heat source tracer test adopts the form of single-point delivery and multi-point monitoring. The hot water source delivery should be implemented in the monitoring point of the center borehole of the face. The temperature of the heat source used should be 40-60 degrees Celsius higher than the ambient background temperature. After the injection time is half of the total monitoring time, the temperature of the monitoring port closest to the injection port should rise by more than 20 degrees Celsius. After the injection of the heat source tracer test is stopped, the subsequent test can be carried out only after the temperature of the monitoring hole of the previous test is guaranteed to return to the initial background temperature. The temperature change of less than 5% is regarded as recovery.

4. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1 is characterized in that: In step (2), the hydraulic test adopts single-point or multi-point pumping and injection, and multi-point continuous monitoring is carried out at the same time, and the experiment is carried out for no less than three times; when the number of monitoring points is less than ten, the pumping and injection experiment is carried out for no less than five times; when the number of monitoring points is more than ten, the number of tests is appropriately reduced; a short monitoring time interval is adopted in the early stage of the test, and the specific interval is shorter than 1% of the total time of the pumping and injection test. As much early head data as possible is collected, and the test is stopped after the head reaches stability, and a quarter of the total time of the designed test is selected as the early stage of pumping and injection.

5. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1 is characterized in that: In step (3), the curve smoothing and noise elimination filtering process is carried out to extract the tracer travel time and head drop data, which specifically includes the following steps: (31) Research and design a smooth low-pass filter function with window weights that changes with variance, carry out pre-denoising, and further use the adaptive Kalman filter method based on conventional Kalman filtering and strong tracking Kalman filtering algorithms to adjust the covariance matrix and gain matrix of the state prediction error in real time according to the dynamic changes of noise, and carry out deep filtering and denoising; (32) For smooth low-pass filtering, a Gaussian filter function is used to denoise the temperature tracer penetration curve. The Gaussian function variance is determined based on the relative size of the data variance in the window and the total variance of the entire data set. When the data variance in the window is greater than the total variance of the data set, the Gaussian function variance is selected as 1.5, otherwise it is selected as 0.

5. Then, combined with the Gaussian function variance, the weight vector of each data in the window is determined, and the pre-denoising temperature curve result is obtained based on the weight and the measured temperature data; (33) For smooth low-pass filtering, a fifth-order or sixth-order polynomial function is used to fit and denoise the water head curve, and the order of the polynomial used must be smaller than the window length. Then, the least squares fitting is performed on the given high-order polynomial, and the pre-denoised water head curve is obtained by stepwise moving average. (34) Assuming that the system noise and the observation noise are uncorrelated, the system noise variance and the measurement noise variance are estimated and corrected by combining the filtering itself, and the time-varying attenuation factor l is introduced. k+1 According to the dynamic changes of noise, the prediction covariance matrix U in the conventional Kalman filter is modified k+1,k , adjust the covariance matrix and gain matrix of the state prediction error in real time; Among them, Trace[] is to find the matrix trace, I k is the identity matrix, is the state transfer matrix, U k is the filter covariance matrix, V k-1 is the residual covariance matrix, Z k is the new information sequence, R k is the system noise covariance matrix, Q k is the measurement noise covariance matrix; (35) Monte Carlo analysis was performed on the denoised temperature tracer penetration curve to extract the peak travel time data. Drawdown data were extracted from the denoised hydraulic head curve. One quarter of the total time after the start of the test was considered the early stage of the pumping test, and the time from the cessation of pumping to the complete recovery of the hydraulic head was considered the late stage of the pumping test. The rest of the time was considered the middle stage. Combining the early, middle and late stage data, no less than five drawdown data were selected from each group of hydraulic head curves, with at least two data selected from each of the early and late stages.

6. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1, characterized in that: In step (4), determine the peak time t of the temperature function derivative B With peak travel time t D The quantitative relationship between the delay factor and t B =δt D (5) In the above formula, and are the first-order and second-order time derivatives of the heat source temperature signal, t B and t D are the derivative peak time and peak time respectively, δ is the conversion factor, C is the convective heat transfer coefficient, D L is the thermal diffusivity; For the travel time dataset after correction, Monte Carlo analysis was carried out to determine different travel time means and standard deviations, including t-5%, t-25%, t-50%, t-75% and t-100% as travel time datasets, and iterations were carried out in turn based on multiple sets of datasets.

7. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1 is characterized in that: In step (5), the heat source temperature tracer migration process is regarded as a series of solute concentration wavefront expansion forms. It is assumed that the initial slowness parameter field is homogeneously distributed. According to the scale of the research domain and the inversion accuracy, the model domain is established and meshed. The medium characteristics are reconstructed iteratively based on the travel time data set and the joint iterative reconstruction algorithm. The slowness correction increment of each unit is determined by combining the ray path lengths of different units. All units are weighted corrected one by one to match the measured ray projection. It is iterated until the travel time difference is less than the predetermined threshold. Then the slowness field is converted into the permeability parameter distribution, which is specifically converted into: Where R is the thermal resistivity, φ is the porosity, i is the local hydraulic gradient, and K is the permeability coefficient. Combined with the multiple travel time data sets obtained in step (6), a parameter sensitivity analysis is performed with the same number of iterations. By comparing the errors between the observed values ​​and the reconstructed values, the parameter distribution obtained by the best reconstructed travel time is selected.

8. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1 is characterized in that: In step (5), when the number of iterations is the same, a parameter sensitivity analysis is performed to select the parameter distribution obtained by the best reconstruction travel time by comparing the errors between the observed values ​​and the reconstructed values.

9. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic combined tomographic inversion according to claim 1, characterized in that: In step (6), the statistical parameters of the hydrogeological parameters, including the mean and variance, are estimated based on the parameter distribution obtained by iterative tracer travel time. Subsequently, the parameters of the model domain are set to a non-homogeneous field that satisfies the Gaussian distribution, or the parameter distribution obtained by iterative tracer travel time is used as the initial field. In combination with coring and geophysical exploration, the statistical characteristics of the parameter field are estimated. At the same time, the relevant scale is determined based on the stratigraphic conditions obtained from the coring.

10. The method for advanced prediction of disaster-causing structures based on heat source tracing and hydraulic joint tomography inversion according to claim 1, characterized in that: In step (7), based on the simultaneous continuous linear estimation algorithm, the geostatistical inversion analysis of the head data is performed, and the conditional effective head field is solved in combination with the initial parameter distribution to determine the mean and variance of the hydrogeological parameters, as well as the unconditional covariance and cross covariance of the parameters and the head response, and the relevant contribution coefficient information; Subsequently, the first hydrogeological parameter estimation is determined, the parameter perturbation of the first iteration and the difference between the observed and simulated hydraulic heads of each group of pumping are obtained, and the conditional means of the parameters and their covariances are updated; Based on the parameter distribution after the first iteration, the conditional effective head field is updated, the head disturbance and related covariance are updated, and the contribution coefficient of each parameter at each unit position under different pumping processes at different times is further calculated and determined. The previous parameter estimation is improved by repeated iterations until the updated parameter conditional covariance or head error meets the threshold setting.

Images