CN114036202A - Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion - Google Patents

Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion Download PDF

Info

Publication number
CN114036202A
CN114036202A CN202111174530.4A CN202111174530A CN114036202A CN 114036202 A CN114036202 A CN 114036202A CN 202111174530 A CN202111174530 A CN 202111174530A CN 114036202 A CN114036202 A CN 114036202A
Authority
CN
China
Prior art keywords
monitoring
test
heat source
disaster
water
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111174530.4A
Other languages
Chinese (zh)
Inventor
王欣桐
许振浩
朱其志
张瑨
林鹏
潘东东
王朝阳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202111174530.4A priority Critical patent/CN114036202A/en
Publication of CN114036202A publication Critical patent/CN114036202A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24Querying
    • G06F16/245Query processing
    • G06F16/2458Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2462Approximate or statistical queries
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29Geographical information databases
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/08Construction

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computational Linguistics (AREA)
  • General Health & Medical Sciences (AREA)
  • Software Systems (AREA)
  • Fuzzy Systems (AREA)
  • Remote Sensing (AREA)
  • Health & Medical Sciences (AREA)
  • Economics (AREA)
  • Mathematical Physics (AREA)
  • Human Resources & Organizations (AREA)
  • Marketing (AREA)
  • Primary Health Care (AREA)
  • Strategic Management (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention discloses a disaster-causing structure advanced forecasting method based on heat source tracing and hydraulic united chromatography inversion, which adopts a test monitoring arrangement scheme which is closely combined with advanced geological drilling conditions, is related to tunnel (cave) and tunnel bore diameters, is implemented based on advanced geological drilling performed in underground engineering, has universality for various underground engineering, adopts a heat source tracing method which has good application prospect due to convenient monitoring and no pollution, realizes dynamic adjustment of a covariance matrix and a gain matrix of a state forecasting error based on a filtering mode combining smooth low-pass filtering and self-adaptive Kalman filtering, has the capability of tracking a mutation state, effectively avoids filtering divergence, fully utilizes data information obtained by the traditional hydrogeological test, develops hydrogeological united chromatography, can obviously reduce the multiple resolvability of inverse problems, and realizes interpretation of contained disaster-causing structure water-rich non-redundant information, provides theoretical basis for accurately detecting the disaster-causing structure in front of the project in advance.

Description

Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion
Technical Field
The invention relates to the technical field of advanced forecasting of disaster-causing structures in the field of underground engineering, in particular to an advanced forecasting method of disaster-causing structures based on heat source tracing and hydraulic combined chromatography inversion.
Background
With the shift of the construction emphasis of the hydraulic and hydroelectric engineering and the railway and highway engineering to the western mountainous areas under the complex terrain and geological conditions, the underground engineering construction has the characteristics of large burial depth, complex geology and frequent disasters, and is easy to encounter geological disasters such as water inrush, mud inrush and the like in the engineering construction. The poor geologic body developed in rock mass has very obvious disaster-causing effect, and comprises karst cave, karst pipeline, underground river and other karst poor geologic bodies, natural and artificial caves, fault fracture zone, strong weathering zone and other special rock-soil bodies. Under the combined action of water head pressure and construction disturbance, the water-resisting rock stratum or the water-blocking structure may be unstable, which further induces water and mud bursting disasters, and becomes a main cause of the water and mud bursting disasters. Underground engineering construction needs to accurately find out a water-rich disaster-causing structure existing in front of an engineering in advance to prevent the disaster.
The traditional exploration means applied to the field of geophysical is one of means commonly adopted for carrying out disaster source detection, the obtained parameters such as conductivity and permeability are complex and fuzzy in the conversion process with the formation hydraulic characteristics, and meanwhile, related conversion is often related to specific formation conditions, so that the universal reference is difficult to provide for specific engineering. The method is combined with the latest research in the field of hydrogeology, data information obtained by traditional hydrogeology tests such as a tracing test and a water pumping and injecting test which are carried out on underground engineering contains a plurality of pieces of non-redundant information describing the water richness of the disaster-causing structure, and combined inversion is carried out by combining water head depthkeeping and tracing data, so that a new idea is provided for the advanced detection research of the disaster-causing structure. In recent years, temperature is paid attention as a novel tracer which is low in cost, convenient and pollution-free, the rapid development of a double-seal layered well formation technology and a distributed optical fiber temperature sensor also creates possibility for in-situ real-time measurement with high precision and low cost, and technical support is provided for monitoring water head and temperature tracing data. Under the background, a joint inversion advanced prediction method suitable for disaster-causing structures is urgently needed to be researched by combining heat source tracing and hydraulic information.
Disclosure of Invention
The invention aims to solve the technical problem of providing an advanced disaster-causing structure forecasting method based on heat source tracing and hydraulic combined chromatography inversion, which is used for acquiring hydrogeological parameter distribution and exploring a water-rich disaster-causing structure position and providing guarantee for engineering safety construction.
In order to solve the technical problem, the invention provides an advanced disaster-causing structure forecasting method based on heat source tracing and hydraulic combined chromatography inversion, which comprises the following steps:
(1) based on geological data, geophysical prospecting and drilling data of a specific engineering research area, aiming at a water-rich structure fracture zone, a karst development section and a geophysical abnormal area in short-medium distance forecasting, the length, depth and inclination angle of a pilot borehole performed by advanced geological detection are combined, and a monitoring area in front of a tunnel face is defined; secondly, determining a heat source tracing test and layered source supplementing hydraulic test monitoring scheme by combining the advanced drilling arrangement position and quantity of the underground engineering and the tunnel and roadway hole diameters, wherein the monitoring scheme comprises the monitoring group number, the quantity of each group and the monitoring interval;
(2) according to the determined monitoring area and the test monitoring scheme, a heat source tracing test and a water pumping and injecting hydraulic test are carried out, and a tracing penetration curve and a water head evolution curve are monitored and recorded;
(3) carrying out curve smoothing denoising filtering processing, and extracting tracing travel time and water head depth reduction data;
(4) considering that the thermal diffusion effect and the dispersion effect exist in the temperature propagation process, a delay factor is determined based on the fastest change moment of the first derivative of the temperature penetration time curve, and deviation correction during tracing travel is achieved. In order to avoid the uncertainty of selection in the early travel, the diagnosis and analysis of the early travel are further carried out aiming at the corrected travel time data set;
(5) carrying out iteration to realize medium characteristic reconstruction based on the travel time data set and a joint iteration reconstruction algorithm, and selecting the optimal reconstruction travel time and estimation parameter distribution;
(6) and determining initial estimation mean, variance, relevant scales and the like of the aquifer parameters by combining the estimated hydrogeological parameter distribution obtained based on travel as prior information for further carrying out the hydraulic chromatography analysis.
(7) And carrying out hydro-chromatographic inversion analysis by combining a simultaneous continuous linear estimation algorithm, carrying out successive iteration until threshold setting is met, then calculating and determining hydrogeological parameter distribution such as permeability coefficient, water storage coefficient, pressure conduction coefficient and the like, and determining the position and distribution of the water-rich disaster-causing structure by combining a high-permeability zone.
Preferably, in the step (1), a depth range of 1.2-1.5 times of the longest drilled hole in front of the tunnel face is taken as a monitoring area, the depth range is 30-100 m, if the engineering required detection precision is high, the development degree of karst in a research area is high, the water-rich property is strong, the adopted experimental monitoring interval is 0.5-1 times of the hole diameter, and the maximum 2-3 times of the hole diameter is taken as the monitoring interval; if the advanced drilling holes reveal a potential water-rich disaster-causing structure, monitoring points are encrypted as much as possible by combining a plurality of drilling holes within the range of 3 times of the hole diameter in the potential disaster-causing structure area, the monitoring distance is reduced to be 0.3-0.5 times of the minimum distance of the drilling holes, and fine detection is conducted again through testing to determine the scale and the characteristics of the disaster-causing structure.
Preferably, in the step (2), the heat source tracing test adopts a single-point putting and multi-point monitoring mode, the hot water source putting is preferably carried out in a monitoring point drilled in the center of the tunnel face, the temperature of the used heat source is 40-60 ℃ higher than the ambient background temperature, and after the injection duration is half of the total monitoring duration, the temperature of the monitoring port closest to the injection port is preferably raised by more than 20 ℃; after the heat source tracing test stops injecting, after the temperature of the preorder test monitoring hole is ensured to be recovered to the initial background temperature, a subsequent test can be carried out, and the recovery is considered when the temperature change is less than 5%.
Preferably, in the step (2), the hydraulic test adopts single-point or multi-point pumping injection, meanwhile, multi-point continuous monitoring is carried out, and the test is carried out for not less than three times; when the number of monitoring points is less than ten, the water pumping and injecting experiment is carried out for not less than five times; when the number of monitoring points is more than ten, the test times are properly reduced; short monitoring time intervals are adopted in the early stage of the test, the specific intervals are shorter than 1% of the total time of the water pumping and injecting test, early water head data are collected as much as possible, the test is stopped after the water head is stabilized, and one fourth of the total time of the design test is selected as the early stage of water pumping and injecting.
Preferably, in the step (3), the curve smoothing denoising filtering processing is carried out, and the step of extracting the tracing travel time and water head depth reduction data specifically comprises the following steps:
(31) a smooth low-pass filtering function with window weighting changing along with variance is researched and designed, pre-denoising processing is carried out, a self-adaptive Kalman filtering method based on a conventional Kalman filtering algorithm and a strong tracking Kalman filtering algorithm is further carried out, a covariance matrix and a gain matrix of a state prediction error are adjusted in real time according to dynamic change of noise, and deep filtering denoising is carried out;
(32) for smooth low-pass filtering, denoising the temperature tracing penetration curve by adopting a Gaussian filter function, determining the variance of the Gaussian function according to the relative size of the data variance in the window and the total variance of the whole data set, selecting the variance of the Gaussian function to be 1.5 when the data variance in the window is larger than the total variance of the data set, and otherwise, selecting the variance of the Gaussian function to be 0.5, further determining the weight vector of each data in the window by combining the variance of the Gaussian function, and acquiring a denoising temperature curve result according to the weight and the measured temperature data;
(33) aiming at the smooth low-pass filtering, fitting and denoising the waterhead curve by adopting a fifth-order or sixth-order polynomial function, wherein the order of the polynomial is required to be smaller than the length of a window, and then carrying out least square fitting on a given high-order polynomial, and gradually moving and averaging to obtain a pre-denoising waterhead curve;
(34) assuming that system noise and observation noise are not correlated with each other, estimating and correcting system noise variance and measured noise variance by combining filtering, and introducing time-varying attenuation factor lk+1Based on dynamic variation of noise, correcting in conventional Kalman filteringPrediction covariance matrix Uk+1,kAdjusting a covariance matrix and a gain matrix of the state prediction error in real time;
Figure BDA0003294487510000031
Figure BDA0003294487510000033
wherein, Trace [ 2 ]]To solve the matrix trace, IkIs a matrix of the units,
Figure BDA0003294487510000032
being a state transition matrix, UkFor filtering covariance matrix, Vk-1As a residual covariance matrix, ZkAs a neoformation sequence, RkIs a system noise covariance matrix, QkMeasuring a noise covariance matrix;
(35) carrying out Monte Carlo analysis to extract peak travel time data according to the temperature tracing penetration curve after noise elimination; extracting depth reduction data aiming at the water head curve after noise elimination, taking a quarter of total duration after starting the test as the early stage of the water pumping test, taking the water pumping stop till the water head is completely recovered as the late stage of the water pumping test, and taking the rest time as the middle stage; and combining the early, middle and late stage data, selecting not less than five depth reduction data in each group of water head curve, wherein at least two data are respectively selected in the early and late stage time.
Preferably, in step (4), the time t at which the derivative of the temperature function peaks is determinedBAnd peak travel time tDFurther obtaining said conversion factor as
Figure BDA0003294487510000041
Figure BDA0003294487510000042
Figure BDA0003294487510000043
tB=δtD (5)
In the above formula, the first and second carbon atoms are,
Figure BDA0003294487510000044
and
Figure BDA0003294487510000045
first and second time derivatives, t, of the heat source temperature signalBAnd tDRespectively, the derivative peak time and the peak time, delta is a conversion factor, C is a convective heat transfer coefficient, and DLIs the thermal diffusivity;
and (3) aiming at the corrected travel time data set, carrying out Monte Carlo analysis to determine specific different travel time mean values and standard deviations thereof, wherein t-5%, t-25%, t-50%, t-75% and t-100% are used as travel time data sets, and carrying out iteration in sequence according to a plurality of groups of data sets.
Preferably, in the step (5), the heat source temperature tracing and transporting process is regarded as a series of wave front expansion forms of solute concentration, an initial slowness parameter field is assumed to be in homogeneous distribution, and a model domain is established and is subjected to grid division according to the research domain scale and the inversion precision; iteration is carried out in sequence to reconstruct medium characteristics based on a travel time data set and a combined iterative reconstruction algorithm, slowness correction increments of all units are determined by combining the path lengths of rays of different units, successive weighted correction is carried out on all the units to match measured ray projections, successive iteration is carried out until the travel time difference is smaller than a preset threshold value, and then a slowness field is converted into penetration parameter distribution, specifically:
Figure BDA0003294487510000046
wherein R is thermal resistivity, phi is porosity, i is local hydraulic gradient, and K is permeability coefficient; and (4) combining the multiple groups of travel time data sets obtained in the step (6), carrying out parameter sensitivity analysis under the condition of the same iteration times, and selecting the optimal parameter distribution obtained by reconstructing travel time by comparing the errors of the observed value and the reconstructed value.
Preferably, in step (5), under the condition that the iteration times are the same, parameter sensitivity analysis is carried out, and the optimal parameter distribution obtained during reconstruction travel is selected by comparing errors of the observed value and the reconstructed value.
Preferably, in the step (6), the statistical parameters of the hydrogeological parameters are estimated according to the parameter distribution obtained by iteration during tracing travel, including the mean value and the variance, and then the parameters of the model domain are all set to be a heterogeneous field meeting Gaussian distribution, or the parameter distribution obtained by iteration during tracing travel is used as an initial field, and the statistical characteristics of the parameter field are estimated by combining a drilling coring and geophysical prospecting means; meanwhile, the relevant scale is determined by combining the stratum condition obtained by the drilling coring condition.
Preferably, in the step (7), performing geostatistical inversion analysis of the water head data based on a simultaneous continuous linear estimation algorithm, solving a conditional effective water head field by combining initial parameter distribution, and determining a hydrogeological parameter mean value and variance, and unconditional covariance and cross covariance of parameters and water head response thereof, and related contribution coefficient information;
then, determining first hydrogeological parameter estimation, obtaining parameter disturbance of first iteration and differences between observation water heads and simulated water heads of all groups of pumped water, and updating parameter condition mean values and covariance thereof;
updating the effective water head field under the condition based on the parameter distribution after the first iteration, updating the water head disturbance and the related covariance, further calculating and determining different water pumping processes of the water head disturbance at different moments, and improving the previous parameter estimation by repeating the iteration until the updated parameter condition covariance or water head error meets the threshold setting.
The invention has the beneficial effects that: the implementation basis of the method is that advanced geological drilling is performed in underground engineering, repeated hole distribution operation is not needed, and the method has universality on the underground engineering under different geological conditions; meanwhile, the adopted test monitoring arrangement scheme is closely combined with the drilling hole length, depth, inclination angle, position and quantity conditions of advanced geological drilling, is related to the tunnel (cave) and roadway hole diameter, and can dynamically adjust the specific monitoring scheme according to the geological structure positions and distribution conditions of the front water guide fault, water-filled cavern and the like disclosed by drilling; the adopted heat source tracing method has good application prospect due to convenient monitoring and no pollution, and adopts a filtering mode combining smooth low-pass filtering and self-adaptive Kalman filtering, so that the dynamic adjustment of a covariance matrix and a gain matrix of a state prediction error can be realized, the sudden change state tracking capability is realized, and filtering divergence is effectively avoided; by developing the joint chromatographic inversion of the hydrogeological data, the conventional data information obtained by utilizing hydrogeological tests such as a tracing test, a water pumping and injecting test and the like is fully fused, the multi-resolvability of the inverse problem is obviously reduced, the non-redundant information of the water richness of the disaster-causing structure contained in the multi-source hydrogeological data is interpreted, and a theoretical basis is provided for accurately detecting the disaster-causing structure existing in front of the engineering in advance.
Drawings
FIG. 1 is an overall schematic diagram of a heat source tracing and hydraulic test monitoring arrangement for underground engineering.
FIG. 2 is a schematic exploded view of a subsurface engineering heat source tracing and hydraulic test monitoring arrangement of the present invention.
Detailed Description
A disaster-causing structure advanced forecasting method based on heat source tracing and hydraulic combined chromatography inversion comprises the following steps:
(1) based on geological data, geophysical prospecting and drilling data of a specific engineering research area, aiming at a water-rich structure fracture zone, a karst development section and a geophysical abnormal area in short-medium distance forecasting, the conditions of the length, depth and inclination angle of a lead borehole performed by advanced geological detection are combined, and a depth range which is 1.2-1.5 times of the longest borehole in front of a tunnel face is defined as a monitoring area, and is usually in a range of 30-100 meters. And then, determining a monitoring scheme of a heat source tracing test and a layered source supplementing hydraulic test by combining the arrangement position and the number of the advanced drilling of the underground engineering, the tunnel diameter and the tunnel diameter, and the like, wherein the monitoring scheme comprises the number of monitoring groups, the number of each group and the monitoring interval.
In the underground engineering such as tunnel, tunnel and tunnel, combine the construction design of advance drilling, arrange the monitoring point in the three-dimensional space of engineering tunnel face the place ahead, figure 1 and figure 2 are one of the example schemes, can combine the quantity of advance drilling to confirm the monitoring group number, confirm the concrete monitoring point number of every group according to the hole diameter size of tunnel, tunnel and tunnel, etc., if the required detection precision of engineering is higher and research district karst development degree is high, the water-rich is strong, the experimental monitoring interval that can adopt is 0.5-1 times the hole diameter, otherwise can adopt 2-3 times the hole diameter as the monitoring interval at most. In addition, if the advanced drilling holes reveal a potential water-rich disaster-causing structure, monitoring points are encrypted as much as possible by combining a plurality of drilling holes within 3 times of the hole diameter range near the potential disaster-causing structure area, the monitoring distance can be considered to adopt 0.3-0.5 times of the minimum distance of the drilling holes, and fine detection is carried out by re-testing to determine the scale and the characteristics of the disaster-causing structure;
(2) and according to the determined monitoring area and the test monitoring scheme, performing a heat source tracing test and a water pumping and injecting hydraulic test, and monitoring and recording a tracing penetration curve and a water head curve.
The heat source tracing test can adopt a single-point feeding and multi-point monitoring mode, the feeding of the hot water source is preferably selected in a monitoring point of a central drilling hole of a tunnel face, the temperature of the used heat source is 40-60 ℃ higher than the ambient background temperature, and after the injection duration is half of the total monitoring duration, the temperature of a monitoring port closest to the injection port is preferably raised by more than 20 ℃. After the injection is stopped, after the temperature of the monitoring hole is recovered to the initial background temperature after the preorder test is ensured, the subsequent test can be carried out, and generally, the temperature change is considered to be recovered when the temperature change is less than 5 percent.
The hydraulic test adopts single-point or multi-point pumping injection, simultaneously multi-point continuous monitoring is carried out, and when the number of monitoring points is less than ten, the pumping injection test is carried out for not less than five times; when the number of monitoring points is more than ten, the number of tests can be reduced appropriately. Because the early stage of the test is fast in water head change, short monitoring time intervals are adopted in the early stage of the test, the specific intervals are shorter than 1% of the total time of the water pumping and injecting test, the early stage water head data are collected as much as possible, the test is stopped after the water head is stabilized, and the time length of one fourth of the total time of the design test can be considered to be selected as the early stage of the water pumping and injecting.
(3) And carrying out curve smoothing denoising filtering processing, and extracting data such as tracing travel time, water head depth reduction and the like. A smooth low-pass filtering function with window weighting changing along with variance is researched and designed, pre-denoising processing is carried out, and then depth filtering denoising is carried out by combining with an adaptive Kalman filtering method;
for smooth low-pass filtering, a Gaussian filter function can be considered for denoising a temperature tracing penetration curve, in order to ensure that the characteristics of the curve are kept as much as possible while denoising is ensured, the variance of the Gaussian function can be determined according to the relative size of the data variance in a window and the total variance of the whole data set, when the data variance in the window is larger than the total variance, the variance of the Gaussian function is selected to be 1.5, and otherwise, the variance of the Gaussian function is selected to be 0.5. Determining a weight vector of each data in the window by combining the variance of the Gaussian function, and acquiring a pre-denoising result according to the weight and the actually measured temperature data; on the other hand, fitting noise elimination can be carried out on the waterhead curve by adopting a fifth-order or sixth-order polynomial function, the order of the polynomial is required to be smaller than the length of the window, a weighting coefficient is determined according to the ratio of the data variance in the window to the total variance of the whole data set, then least square fitting is carried out on the given high-order polynomial, and the waterhead curve is gradually moved and averaged to obtain the pre-noise elimination waterhead curve.
Aiming at the self-adaptive Kalman filtering, assuming that system noise and observation noise are not correlated with each other, estimating and correcting system noise variance and measurement noise variance are carried out by combining the filtering, and a time-varying attenuation factor l is introducedk+1Correcting a prediction covariance matrix U in conventional Kalman filtering according to dynamic changes of noisek+1,kAnd adjusting the covariance matrix and the 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 change of external noise, so that the filter has the capability of tracking the abrupt change state and avoids filtering divergence.
Figure BDA0003294487510000071
Figure BDA0003294487510000072
UkIn order to filter the covariance matrix,
Figure BDA0003294487510000073
being a state transition matrix, QkFor measuring the covariance matrix of noise, I is the identity matrix, Trace [ ]]To solve the trace of the matrix, ZkAs a neoformation sequence, RkIs a system noise covariance matrix, Vk-1Is a residual covariance matrix.
And (4) carrying out Monte Carlo analysis to extract peak travel time data according to the temperature tracing penetration curve after noise elimination. And (3) extracting the depth reduction data aiming at the water head curve after noise elimination, generally regarding a quarter of total duration after the test is started as the early stage of the water pumping test, regarding the time after the water pumping is stopped until the water head is completely recovered as the late stage of the water pumping test, and regarding the rest time as the middle stage. And combining the early, middle and late stage data, selecting not less than five depth reduction data in each group of water head curve, wherein at least two data are respectively selected in the early and late stage time.
(4) Considering the thermal diffusion effect and the dispersion effect in the temperature propagation process, a delay factor is determined based on the fastest change moment of the first derivative of the temperature penetration time curve, and deviation correction in tracing travel is achieved. In order to avoid the uncertainty of selection in the early travel, diagnosis and analysis in the early travel are further carried out aiming at the corrected travel time data set.
According to the heat convection diffusion control equation, determining the quantitative relation between the temperature function derivative peak time and the peak travel time, and further obtaining the conversion factor of
Figure BDA0003294487510000081
Figure BDA0003294487510000082
Figure BDA0003294487510000083
tB=δtD (5)
In the above formula, the first and second carbon atoms are,
Figure BDA0003294487510000084
and
Figure BDA0003294487510000085
first and second time derivatives, t, of the heat source temperature signalBAnd tDRespectively, the derivative peak time and the peak time, delta is a conversion factor, C is a convective heat transfer coefficient, and DLIs the thermal diffusivity.
And (3) aiming at the corrected travel time data set, carrying out Monte Carlo analysis to determine specific different travel time mean values and standard deviations thereof, wherein t-5%, t-25%, t-50%, t-75% and t-100% are used as travel time data sets, and carrying out iteration in sequence according to a plurality of groups of data sets.
(5) And (3) carrying out iteration to reconstruct the medium characteristics based on the travel time data set and the joint iteration reconstruction algorithm, and selecting the optimal reconstruction travel time and the estimation parameter distribution thereof.
The method specifically comprises the steps of regarding a heat source temperature tracing migration process as a series of wave front expansion forms of solute concentration, assuming that an initial slowness parameter field is in homogeneous distribution, establishing a model domain and performing grid division according to the research domain scale and inversion accuracy.
Then, iteration is carried out in sequence to carry out medium feature reconstruction based on a temperature travel time data set and a combined iteration reconstruction algorithm, slowness correction increments of all units are determined by combining the path lengths of rays of different units, successive weighting correction is carried out on all the units to match measured ray projections, successive iteration is carried out until the travel time difference is smaller than a preset threshold value, then a slowness field is converted into penetration parameter distribution, and specific conversion can refer to:
Figure BDA0003294487510000086
wherein R is thermal resistivity, phi is porosity, i is local hydraulic gradient, and K is permeability coefficient. And under the condition of the same iteration times, carrying out parameter sensitivity analysis, and selecting the optimal parameter distribution obtained during reconstruction travel by comparing the errors of the observed value and the reconstructed value.
(6) And determining information such as initial estimation mean, variance and related scale of aquifer parameters by combining the estimated hydrogeological parameter distribution obtained based on travel.
In some implementation cases, the estimation of hydrogeological parameter statistical parameters including mean value, variance and the like is carried out according to parameter distribution obtained by iteration during tracing travel, and then parameters of a model domain are set to be a heterogeneous field meeting Gaussian distribution; or parameter distribution obtained by iteration during tracing travel can be used as an initial field, and the statistical characteristics of the parameter field are estimated by combining means such as drilling coring and geophysical prospecting. Furthermore, the determination of the relevant dimensions should take into account the formation conditions resulting from the coring of the borehole.
(7) And carrying out hydrogeological parameter inversion analysis by combining with a simultaneous continuous linear estimation algorithm, carrying out successive iteration until threshold setting is met, then calculating and determining parameter distribution such as permeability coefficient, water storage coefficient, pressure conduction coefficient and the like, and determining the position and distribution of the water-rich disaster-causing structure by combining with a high permeability zone.
Specifically, based on a simultaneous continuous linear estimation algorithm, performing geostatistical inversion analysis on water head data, solving a conditional effective water head field by combining initial parameter distribution, and determining information such as a hydrogeological parameter mean value, a hydrogeological parameter variance, unconditional covariance and cross covariance of parameters and water head response, related contribution coefficients and the like.
And then, determining first hydrogeological parameter estimation, obtaining parameter disturbance of first iteration and differences between the observation water heads and the simulated water heads of all groups of pumped water, and updating parameter condition mean values and covariance thereof. Updating the effective water head field under the condition based on the parameter distribution after the first iteration, updating the water head disturbance and the related covariance, further calculating and determining different water pumping processes of the water head disturbance at different moments, and improving the previous parameter estimation by repeating the iteration until the updated parameter condition covariance or water head error meets the threshold setting.

Claims (10)

1. A disaster-causing structure advanced forecasting method based on heat source tracing and hydraulic combined tomography inversion is characterized by comprising the following steps:
(1) based on geological data, geophysical prospecting and drilling data of a specific engineering research area, aiming at a water-rich structure fracture zone, a karst development section and a geophysical abnormal area in short-medium distance forecasting, the length, depth and inclination angle of a pilot borehole performed by advanced geological detection are combined, and a monitoring area in front of a tunnel face is defined; secondly, determining a heat source tracing test and layered source supplementing hydraulic test monitoring scheme by combining the advanced drilling arrangement position and quantity of the underground engineering and the tunnel and roadway hole diameters, wherein the monitoring scheme comprises the monitoring group number, the quantity of each group and the monitoring interval;
(2) according to the determined monitoring area and the test monitoring scheme, a heat source tracing test and a water pumping and injecting hydraulic test are carried out, and a tracing penetration curve and a water head evolution curve are monitored and recorded;
(3) carrying out curve smoothing denoising filtering processing, and extracting tracing travel time and water head depth reduction data;
(4) considering that the thermal diffusion effect and the dispersion effect exist in the temperature propagation process, a delay factor is determined based on the fastest change moment of the first derivative of the temperature penetration time curve, and deviation correction during tracing travel is achieved. In order to avoid the uncertainty of selection in the early travel, the diagnosis and analysis of the early travel are further carried out aiming at the corrected travel time data set;
(5) carrying out iteration to realize medium characteristic reconstruction based on the travel time data set and a joint iteration reconstruction algorithm, and selecting the optimal reconstruction travel time and estimation parameter distribution;
(6) and determining initial estimation mean, variance, relevant scales and the like of the aquifer parameters by combining the estimated hydrogeological parameter distribution obtained based on travel as prior information for further carrying out the hydraulic chromatography analysis.
(7) And carrying out hydro-chromatographic inversion analysis by combining a simultaneous continuous linear estimation algorithm, carrying out successive iteration until threshold setting is met, then calculating and determining hydrogeological parameter distribution such as permeability coefficient, water storage coefficient, pressure conduction coefficient and the like, and determining the position and distribution of the water-rich disaster-causing structure by combining a high-permeability zone.
2. The disaster-causing structure advanced forecasting method based on heat source tracing and hydraulic power combined tomography inversion as claimed in claim 1, characterized in that in step (1), the depth range of 1.2-1.5 times of the longest borehole in front of the tunnel face is taken as a monitoring area, the depth range is 30 m to 100 m, if the detection precision required by engineering is higher, the development degree of karst in a research area is high, the water-rich property is strong, the adopted test monitoring interval is 0.5-1 times of the hole diameter, and 2-3 times of the hole diameter is maximally adopted as the monitoring interval; if the advanced drilling holes reveal a potential water-rich disaster-causing structure, monitoring points are encrypted as much as possible by combining a plurality of drilling holes within the range of 3 times of the hole diameter in the potential disaster-causing structure area, the monitoring distance is reduced to be 0.3-0.5 times of the minimum distance of the drilling holes, and fine detection is conducted again through testing to determine the scale and the characteristics of the disaster-causing structure.
3. The disaster-causing structure advanced forecasting method based on heat source tracing and hydraulic power combined tomography inversion as claimed in claim 1, characterized in that in step (2), the heat source tracing test adopts a single-point putting and multi-point monitoring mode, the putting of the hot water source is preferably implemented in a monitoring point of a central borehole of a tunnel face, the temperature of the used heat source is 40-60 ℃ higher than the ambient background temperature, and after the injection duration is half of the total monitoring duration, the temperature of the monitoring port nearest to the injection port is preferably raised by more than 20 ℃; after the heat source tracing test stops injecting, after the temperature of the preorder test monitoring hole is ensured to be recovered to the initial background temperature, a subsequent test can be carried out, and the recovery is considered when the temperature change is less than 5%.
4. The disaster-causing structure advanced forecasting method based on heat source tracing and hydraulic combined tomography inversion as claimed in claim 1, characterized in that in step (2), single-point or multi-point pumping injection is adopted in the hydraulic test, multi-point continuous monitoring is carried out simultaneously, and the experiment is carried out not less than three times; when the number of monitoring points is less than ten, the water pumping and injecting experiment is carried out for not less than five times; when the number of monitoring points is more than ten, the test times are properly reduced; short monitoring time intervals are adopted in the early stage of the test, the specific intervals are shorter than 1% of the total time of the water pumping and injecting test, early water head data are collected as much as possible, the test is stopped after the water head is stabilized, and one fourth of the total time of the design test is selected as the early stage of water pumping and injecting.
5. The disaster-causing structure advanced prediction method based on heat source tracing and hydraulic power combined tomography inversion as claimed in claim 1, wherein in the step (3), curve smoothing noise-eliminating filtering processing is carried out, and the steps of extracting tracing travel time and water head depth reduction data specifically comprise the following steps:
(31) a smooth low-pass filtering function with window weighting changing along with variance is researched and designed, pre-denoising processing is carried out, a self-adaptive Kalman filtering method based on a conventional Kalman filtering algorithm and a strong tracking Kalman filtering algorithm is further carried out, a covariance matrix and a gain matrix of a state prediction error are adjusted in real time according to dynamic change of noise, and deep filtering denoising is carried out;
(32) for smooth low-pass filtering, denoising the temperature tracing penetration curve by adopting a Gaussian filter function, determining the variance of the Gaussian function according to the relative size of the data variance in the window and the total variance of the whole data set, selecting the variance of the Gaussian function to be 1.5 when the data variance in the window is larger than the total variance of the data set, and otherwise, selecting the variance of the Gaussian function to be 0.5, further determining the weight vector of each data in the window by combining the variance of the Gaussian function, and acquiring a denoising temperature curve result according to the weight and the measured temperature data;
(33) aiming at the smooth low-pass filtering, fitting and denoising the waterhead curve by adopting a fifth-order or sixth-order polynomial function, wherein the order of the polynomial is required to be smaller than the length of a window, and then carrying out least square fitting on a given high-order polynomial, and gradually moving and averaging to obtain a pre-denoising waterhead curve;
(34) assuming that system noise and observation noise are not correlated with each other, estimating and correcting system noise variance and measured noise variance by combining filtering, and introducing time-varying attenuation factor lk+1Correcting a prediction covariance matrix U in conventional Kalman filtering according to dynamic changes of noisek+1,kAdjusting covariance matrix and gain matrix of state prediction error in real time;
Figure FDA0003294487500000021
Figure FDA0003294487500000031
Wherein, Trace [ 2 ]]To solve the matrix trace, IkIs a matrix of the units,
Figure FDA0003294487500000032
being a state transition matrix, UkFor filtering covariance matrix, Vk-1As a residual covariance matrix, ZkAs a neoformation sequence, RkIs a system noise covariance matrix, QkMeasuring a noise covariance matrix;
(35) carrying out Monte Carlo analysis to extract peak travel time data according to the temperature tracing penetration curve after noise elimination; extracting depth reduction data aiming at the water head curve after noise elimination, taking a quarter of total duration after starting the test as the early stage of the water pumping test, taking the water pumping stop till the water head is completely recovered as the late stage of the water pumping test, and taking the rest time as the middle stage; and combining the early, middle and late stage data, selecting not less than five depth reduction data in each group of water head curve, wherein at least two data are respectively selected in the early and late stage time.
6. The method for forecasting the disaster-causing structure in advance based on the heat source tracing and the hydraulic combined tomography inversion as claimed in claim 1, wherein in the step (4), the peak time t of the derivative of the temperature function is determinedBAnd peak travel time tDFurther obtaining said conversion factor as
Figure FDA0003294487500000033
Figure FDA0003294487500000034
Figure FDA0003294487500000035
tB=δtD (5)
In the above formula, the first and second carbon atoms are,
Figure FDA0003294487500000036
and
Figure FDA0003294487500000037
first and second time derivatives, t, of the heat source temperature signalBAnd tDRespectively, the derivative peak time and the peak time, delta is a conversion factor, C is a convective heat transfer coefficient, and DLIs the thermal diffusivity;
and (3) aiming at the corrected travel time data set, carrying out Monte Carlo analysis to determine specific different travel time mean values and standard deviations thereof, wherein t-5%, t-25%, t-50%, t-75% and t-100% are used as travel time data sets, and carrying out iteration in sequence according to a plurality of groups of data sets.
7. The disaster-causing structure advanced prediction method based on heat source tracing and hydraulic power combined tomography inversion as claimed in claim 1, wherein in step (5), the heat source temperature tracing migration process is regarded as a series of wavefront expansion forms of solute concentration, an initial slowness parameter field is assumed to be in homogeneous distribution, and a model domain is established and grid division is performed according to the scale of a research domain and the inversion accuracy; iteration is carried out in sequence to reconstruct medium characteristics based on a travel time data set and a combined iterative reconstruction algorithm, slowness correction increments of all units are determined by combining the path lengths of rays of different units, successive weighted correction is carried out on all the units to match measured ray projections, successive iteration is carried out until the travel time difference is smaller than a preset threshold value, and then a slowness field is converted into penetration parameter distribution, specifically:
Figure FDA0003294487500000041
wherein R is thermal resistivity, phi is porosity, i is local hydraulic gradient, and K is permeability coefficient; and (4) combining the multiple groups of travel time data sets obtained in the step (6), carrying out parameter sensitivity analysis under the condition of the same iteration times, and selecting the optimal parameter distribution obtained by reconstructing travel time by comparing the errors of the observed value and the reconstructed value.
8. The disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined tomography inversion as claimed in claim 1, wherein in step (5), under the condition of the same iteration times, parameter sensitivity analysis is carried out, and by comparing the error between the observed value and the reconstructed value, the optimal parameter distribution obtained during reconstruction travel is selected.
9. The disaster-causing structure advanced forecasting method based on heat source tracing and hydraulic power combined tomography inversion as claimed in claim 1, characterized in that in step (6), the hydrogeological parameter statistical parameters are estimated according to the parameter distribution obtained by iteration during tracing travel, including the mean value and the variance, and then the parameters of the model domain are all set to be a heterogeneous field satisfying the gaussian distribution, or the parameter distribution obtained by iteration during tracing travel is adopted as an initial field, and the statistical properties of the parameter field are estimated by combining the drilling coring and geophysical prospecting means; meanwhile, the relevant scale is determined by combining the stratum condition obtained by the drilling coring condition.
10. The disaster-causing structure advanced prediction method based on heat source tracing and hydraulic power joint tomography inversion as claimed in claim 1, wherein in step (7), geostatistical inversion analysis of water head data is performed based on a simultaneous continuous linear estimation algorithm, a conditional effective water head field is solved by combining initial parameter distribution, and the mean and variance of hydrogeological parameters, unconditional covariance and cross covariance of parameters and water head response, and related contribution coefficient information are determined;
then, determining first hydrogeological parameter estimation, obtaining parameter disturbance of first iteration and differences between observation water heads and simulated water heads of all groups of pumped water, and updating parameter condition mean values and covariance thereof;
updating the effective water head field under the condition based on the parameter distribution after the first iteration, updating the water head disturbance and the related covariance, further calculating and determining different water pumping processes of the water head disturbance at different moments, and improving the previous parameter estimation by repeating the iteration until the updated parameter condition covariance or water head error meets the threshold setting.
CN202111174530.4A 2021-10-09 2021-10-09 Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion Pending CN114036202A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111174530.4A CN114036202A (en) 2021-10-09 2021-10-09 Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111174530.4A CN114036202A (en) 2021-10-09 2021-10-09 Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion

Publications (1)

Publication Number Publication Date
CN114036202A true CN114036202A (en) 2022-02-11

Family

ID=80134760

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111174530.4A Pending CN114036202A (en) 2021-10-09 2021-10-09 Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion

Country Status (1)

Country Link
CN (1) CN114036202A (en)

Similar Documents

Publication Publication Date Title
RU2489735C2 (en) Describing underground structure with function-based iterative inversion
CA2618555C (en) Method and system for pre-drill pore pressure prediction
Zuo et al. Geothermal regime and source rock thermal evolution in the Chagan sag, Inner Mongolia, northern China
Peters et al. Criteria to determine borehole formation temperatures for calibration of basin and petroleum system models
CN111045114B (en) Method for identifying and positioning favorable sand bodies of basalt coverage area sandstone-type uranium deposit mineralization
CN104597490A (en) Multi-wave AVO reservoir elastic parameter inversion method based on precise Zoeppritz equation
Landa et al. Reservoir characterization constrained to well-test data: a field example
CN109799540B (en) Volcanic rock type uranium deposit magnetic susceptibility inversion method based on geological information constraint
CN110737018B (en) Method for modeling anisotropy of VSP seismic data
CN112363210B (en) Coal thickness quantitative prediction method based on transmission groove wave velocity and attenuation coefficient joint inversion
CN112799140B (en) Permeability estimation method based on natural potential inversion
Behm Feasibility of borehole ambient noise interferometry for permanent reservoir monitoring
Desper et al. Accurate water‐table depth estimation using seismic refraction in areas of rapidly varying subsurface conditions
CN114036202A (en) Disaster-causing structure advanced prediction method based on heat source tracing and hydraulic combined chromatography inversion
CN106990433B (en) A kind of recognition methods of the small erosion channel in massif
Dong et al. Diagnosis of concentrated leakage channel embedded in dam base by means of hydraulic tomography
CN112906465B (en) Coal measure stratum acoustic curve reconstruction method and system based on stratum factors
Hua et al. Reservoir-induced seismicity in high seismicity region—a case study of the Xiaowan reservoir in Yunnan province, China
CN109901221B (en) Seismic data anisotropy modeling method based on dynamic correction velocity parameter
CN111965720A (en) Method for acquiring hydraulic conductivity coefficient based on ground-well combination
CN114624779A (en) Pre-stack multi-parameter inversion method for balanced model constraint
CN113495293B (en) Reservoir fluid prediction method and device
CN113267808B (en) Amplitude compensation method and device
Abdrakhimova Improving groundwater flow model parametrization techniques
Zeng et al. Deep-thin volcanic reservoirs characterization using spatial constraint well log-seismic joint inversion

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination