A hydrogeological parameter space-time heterogeneous geological statistics inversion method

Abstract

This invention discloses a geostatistical inversion method for spatiotemporal heterogeneity of hydrogeological parameters. The method includes: acquiring time series of water level, water quality, and pumping flow rate of monitoring wells, as well as the top and bottom depths and interlayer thicknesses of the screen pipe covering the monitoring wells; aligning timestamps to form a hybrid observation dataset; subdividing the layers according to the top and bottom depths of the layers, establishing a layer response sequence using pumping flow rate as an excitation, and matching it with water level and water quality to obtain a contribution weight sequence and a phase difference sequence; calculating an exchange discriminant based on the contribution weight sequence and phase difference sequence and comparing it with a threshold, outputting the depth range of candidate interlayer influence zones if the threshold is exceeded; establishing a vertical inversion model with candidate zones as constraints, synthesizing weighted observations according to contribution weights to participate in fitting, and inverting vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity; changing the pumping flow rate recalculation prediction error and comparing it with the error threshold to identify misaligned fitting.

Description

Technical Field

[0001] This invention relates to the field of data storage technology, and specifically to a method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters. Background Technology

[0002] In multi-layered aquifer systems or systems with weak interlayer connectivity, the screens of monitoring wells often cover multiple layers, resulting in mixed observations of water level and water quality time series, both contributed by multiple layers. Existing inversion processes typically use these mixed observations directly as input models for single-layer observations, lacking a decomposition mechanism to establish response sequences for each subdivided layer based on pumping flow time series and to calculate contribution weight sequences and phase difference sequences. They also lack a constraint mechanism to calculate exchange discriminant quantities based on contribution weight sequences and phase difference sequences, locate candidate segments affecting interlayers, and use these candidate segments to constrain the inversion of vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity. Furthermore, they lack a verification mechanism to recalculate prediction errors under varying pumping flow time series conditions and compare them with preset error thresholds to identify misaligned fits. Consequently, the phase differences and mixed response differences caused by interlayer exchange or wellbore short-circuit are converted into horizontal parameter anomalies, resulting in a valid parameter field fit but invalid layer correspondences, distorting the basis for subsequent risk assessments and engineering layer selection.

[0003] Therefore, proposing a geostatistical inversion method for spatiotemporal heterogeneity of hydrogeological parameters to solve the difficulties of existing technologies is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0004] The purpose of this invention is to provide a geostatistical inversion method for spatiotemporal heterogeneity of hydrogeological parameters to address the shortcomings of the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters, comprising: S101. Obtain the water level time series and water quality time series of the monitoring well before pumping, during pumping, and during the recovery period. Simultaneously obtain the pumping flow time series of the pumping well, as well as the top and bottom depths and interlayer thickness of the screen pipe covering layer of the monitoring well. Align the above data with a unified timestamp to form a hybrid observation dataset. S102. Based on the top and bottom depths of the screen tube cover layer, the vertical range corresponding to the mixed observation dataset is divided into subdivided layers, and the pumping flow time series is used as the excitation to establish a response sequence for each subdivided layer; each response sequence is matched with the water level time series and the water quality time series to obtain the contribution weight sequence and phase difference sequence of each subdivided layer. S103. Calculate the exchange discriminant based on the contribution weight sequence and the phase difference sequence; compare the exchange discriminant with the corresponding preset threshold. When any discriminant exceeds the corresponding preset threshold, output the depth range of the interlayer influence candidate segment. S104. Establish a vertical inversion model with the depth range of the candidate section affected by the interlayer as a constraint. Combine the water level time series and water quality time series according to the contribution weight sequence to form a weighted observation, form a weighted observation sequence set to participate in the fitting, and invert the vertical permeability, interlayer conductivity and well short-circuit conductivity. Under the condition of changing the pumping flow rate time series, recalculate the prediction error and compare the prediction error with the preset error threshold to identify the misaligned fitting.

[0006] Furthermore, methods for dividing the vertical range corresponding to the mixed observation dataset into subdivided segments based on the top and bottom depths of the sieve tube cover layer include: Obtain the top and bottom depths of the sieve tube covering layer to determine the upper and lower boundary depths of the vertical range, and construct the layer boundary sequence based on the top and bottom depths of each layer. Obtain the segment boundary sequence and sort it by depth. Traverse the two adjacent boundary depth values ​​in the segment boundary sequence, take the shallower boundary depth as the top depth of the basic segment and take the deeper boundary depth as the bottom depth of the basic segment to generate the basic segment table, and remove the basic segments that do not intersect with the start and end depths of the sieve tube. The top and bottom depths of the basic layer segments are obtained from the basic layer segment table, and the thickness of the basic layer segment is calculated by subtracting the top depth from the bottom depth. The thickness of the basic layer segment is then compared with a preset thickness threshold. Basic layers whose thickness exceeds the thickness threshold are subdivided into multiple sub-segments according to equal depth steps.

[0007] Furthermore, methods for establishing response sequences for each sub-segment using pumping flow rate time series as excitation include: Obtain the pumping flow rate time series and use a unified timestamp sequence; For each subdivided segment, the center depth and thickness of the segment are calculated based on the top and bottom depths of the segment. The center depth of the segment is mapped to a response delay according to a monotonic mapping rule, and the mapping result is compared with a delay threshold to limit the mapping result. The thickness of the segment is mapped to a response decay according to a monotonic mapping rule, and the mapping result is compared with a decay threshold to limit the mapping result. The response kernel sequence of the subdivided segment is constructed based on the response delay and response decay. The response kernel sequence is convolved with the pumping flow time series to obtain the response sequence of the subdivided segments.

[0008] Furthermore, methods for obtaining subdivided response sequences by performing convolution operations on the response kernel sequence and the pumping flow time series include: Obtain the response kernel sequence associated with the subdivision segment identifier, and obtain the pumping flow time series associated with the same pumping event identifier and aligned with a unified timestamp sequence; The sampling interval of the response kernel sequence is compared with the sampling interval of the unified timestamp sequence. If the sampling intervals are inconsistent, the response kernel sequence is resampled according to the sampling interval of the unified timestamp sequence, while keeping the delay segment at zero and the decay segment gradually decreasing. For each target moment in the unified timestamp sequence, the pumping flow rate value at the corresponding time step is traced back according to the position of the response kernel sequence and multiplied with the corresponding kernel weight. The results are then accumulated to obtain the subdivided segment response value. When the traced back time step is earlier than the start time of the pumping flow rate time series, the multiplication and accumulation are performed with a zero value. The response kernel sequence consists of multiple kernel weights arranged in time step order. The kernel weight is the value of the response kernel sequence at the corresponding position and is used to characterize the contribution ratio of a unit pumping disturbance at different traced back time steps. The response sequence of each subdivided segment is obtained by repeatedly multiplying and accumulating the responses at all target times.

[0009] Furthermore, methods for obtaining the contribution weight sequence and phase difference sequence of each subdivision segment include: Obtain the response sequence and water level time series of each sub-segment, determine the peak time of the response sequence and the peak time of the water level time series of each sub-segment, and determine the core interval by taking the minimum value of the peak time as the starting point and the maximum value of the peak time as the ending point. Based on the core interval, a main response window is generated at both ends by an extended time step. Within the main response window, the response values ​​of each subdivided segment are summed according to a unified timestamp to obtain the segment response sum. The contribution weight is determined based on the comparison between the sum of the segment response and the zero threshold. When the sum of the segment response is not greater than the zero threshold, the contribution weight of each sub-segment at the corresponding timestamp is reset to zero. When the sum of the segment response is greater than the zero threshold, the sub-segment response value at the corresponding timestamp is divided by the sum of the segment response to obtain the sub-segment contribution weight at that timestamp. The contribution weight sequence of each sub-segment is output along the unified timestamp. Based on the main response time window, the peak time of the response sequence of the subdivided segment and the peak time of the water level time sequence are determined respectively. The difference between the two peak times and the sampling interval are converted into the step difference. The phase difference sequence is obtained by repeatedly converting along the time axis by sliding time window. The water quality time series includes conductivity and turbidity sequences. The conductivity and turbidity sequences are obtained and adjacent difference sequences are calculated separately. The difference values ​​are compared with the water quality difference threshold. Timestamps exceeding the water quality difference threshold are marked as water quality change times. At each water quality change time, the phase difference sequence values ​​are read to calculate the absolute value of the time difference between the response arrival time and the water quality change time. The absolute value of the time difference is compared with the time difference threshold to determine the matching segment. At the water quality change time, the matching segment is adjusted upward according to the preset weight, and the non-matching segment is adjusted downward proportionally, while keeping the sum of the contribution weights of each sub-segment equal to one.

[0010] Furthermore, methods for calculating the exchange discriminant based on the contribution weight sequence and the phase difference sequence include: The exchange discrimination parameters include phase difference increment, weighted mutation rate, and water level and water quality lag. Obtain the sampling interval and full time period length of the unified timestamp sequence, generate a time window index according to the sliding step size, and determine each time window by the start timestamp and the end timestamp; Obtain the subdivided layer set and the top and bottom depths of each layer. Also obtain the contribution weight sequence and phase difference sequence associated with each subdivided layer identifier. Obtain the water level time series and water quality time series associated with the same pumping event identifier. Establish an alignment table by associating the time window index with the timestamp, subdivided layer, contribution weight, and phase difference. The subdivided layer set is composed of subdivided layers. Based on the alignment table, adjacent subdivided segments are sorted by the center depth of the segment to form pairs. The phase difference increment is calculated in each time window and correlated with the boundary depth of the two segments. The absolute value of the difference between adjacent timestamps of the contribution weight sequence is calculated and divided by the sampling interval to obtain the weight mutation rate sequence, which is then correlated with the depth range of the subdivided segment. Based on the alignment table, the adjacent difference sequences of the water level time series and the adjacent difference sequences of the water quality time series are calculated in each time window. The cross-correlation value sequence is calculated and the time step number corresponding to the maximum cross-correlation value is taken as the lag between water level and conductivity and the lag between water level and turbidity. These are associated with the time window identifier and stored to form the lag between water level and water quality.

[0011] Furthermore, when any discrimination value exceeds the corresponding preset threshold, the method for outputting the depth range of the candidate segment affected by the interlayer includes: The phase difference increment, weighted mutation rate, and water level and water quality lag are obtained in each time window. The phase difference increment is compared with the phase difference increment threshold, the weighted mutation rate is compared with the weighted mutation rate threshold, and the water level and water quality lag is compared with the lag threshold. When any discriminant exceeds the corresponding threshold, the discriminant type, discriminant value, time window identifier, and associated subdivision segment identifier or adjacent segment boundary depth are recorded to form an exchange overthreshold set. For each record in the exchange threshold set, when the discriminant type is phase difference increment, the top depth and bottom depth of the two segments are merged based on the boundary depth of adjacent subdivision segments to obtain the candidate segment depth range. When the discriminant type is weighted mutation rate, the top depth and bottom depth of the subdivision segment are output based on the subdivision segment identifier to obtain the candidate segment depth range. When the discriminant type is water level and water quality lag, the top depth and bottom depth of the segment with the largest weight in the time window contribution weight sequence are output based on the water level and bottom depth of the segment to obtain the candidate segment depth range.

[0012] Furthermore, methods for establishing vertical inversion models using the depth range of candidate segments affected by interlayers as constraints include: Obtain the depth range of interlayer impact candidate segments associated with the same pumping event identifier, and obtain the top depth and bottom depth of the subdivided segment set. Based on the intersection relationship between the depth range and the subdivided segments, mark the candidate segment segments and non-candidate segment segments. A layered parameter set for the vertical inversion model is established based on candidate and non-candidate segments. The layered parameter set includes vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity. Vertical permeability is taken independently by segment within the candidate segment and by a preset sharing rule within the non-candidate segment. Based on the boundary depth of the candidate segment, the adjacent layer pairs that cross the boundary depth of the candidate segment are determined. The inter-layer connectivity is set independently for the adjacent layer pairs and set according to the preset sharing rules for the other adjacent layer pairs. Set parameter boundary threshold constraints for wellbore short-circuit conductivity and correlate wellbore short-circuit conductivity with the top and bottom depths of the screen pipe cover section.

[0013] Furthermore, methods for synthesizing water level time series and water quality time series into weighted observations according to contribution weights, forming a weighted observation sequence set to participate in fitting, and inverting vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity include: Obtain the contribution weight sequence associated with each sub-segment identifier and obtain the water level time series and water quality time series associated with the same pumping event identifier and aligned with a unified timestamp; For each timestamp, the water level observation, conductivity observation and turbidity observation are weighted and synthesized according to the contribution weight of each sub-segment to form a weighted observation sequence set, and the weighting process keeps the sum of the contribution weights of each segment equal to one. The vertical inversion model is input with the pumping flow time series and outputs a set of prediction sequences aligned with a unified timestamp sequence, including at least the water level prediction sequence, conductivity prediction sequence and turbidity prediction sequence. The difference sequences between the water level prediction sequence and the weighted water level observation sequence, the difference sequences between the conductivity prediction sequence and the weighted conductivity observation sequence, and the difference sequences between the turbidity prediction sequence and the weighted turbidity observation sequence are calculated separately, and the sum of the squares of each difference sequence is combined as the fitting error index. Under the constraint of parameter boundary threshold, the vertical permeability, interlayer conductivity and wellbore short-circuit conductivity are iteratively updated, and the fitting error index is recalculated until the fitting error index no longer decreases or reaches the iteration number threshold. The vertical permeability, interlayer conductivity and wellbore short-circuit conductivity are then output.

[0014] Furthermore, the method of recalculating the prediction error under varying pumping flow rate time series conditions and comparing the prediction error with a preset error threshold to identify misfitting layers includes: The modified pumping flow time series is constructed by obtaining the pumping flow time series corresponding to another pumping event identifier, or by performing amplitude transformation and time shift on the pumping flow time series under the same pumping event identifier to generate the modified pumping flow time series. While keeping the vertical permeability, interlayer conductivity and wellbore short-circuit conductivity unchanged, the modified pumping flow rate time series is input into the vertical inversion model to output the recalculated water level prediction series and the recalculated water quality prediction series. Calculate the difference sequence between the recalculated water level prediction sequence and the water level time series under the corresponding conditions, and calculate the difference sequence between the recalculated water quality prediction sequence and the water quality time series under the corresponding conditions. Combine the sum of squares of the difference sequences to form the prediction error. The prediction error is compared with a preset error threshold. When the prediction error exceeds the preset error threshold, a misalignment fitting indicator is output. The depth range of the candidate segment that triggered the misalignment fitting indicator and the corresponding stratification parameter values ​​are output and associated and stored.

[0015] The technical effects and advantages of the geostatistical inversion method for spatiotemporal heterogeneity of hydrogeological parameters provided by this invention are as follows: By transforming the mixed observations caused by screen pipe crossing layers into computable inputs, water level time series, water quality time series, and pumping flow time series are acquired simultaneously within the same pumping event. These are then combined with the top and bottom depths of the screen pipe covered layer and the thickness of the interlayer to form a mixed observation dataset. Abrupt points in water level / water quality are marked and processed through threshold comparison to reduce the impact of non-mechanistic spikes on the calculation of discriminants. As a result, observations, disturbances, and vertical sources are integrated at the data level, reducing the risk of attribution bias caused by asynchronous pumping input and well response and missing layer information, thus providing a clear data foundation for subsequent stratified analysis.

[0016] Using the pumping flow time series as an excitation, the system is divided into subdivided segments according to the screen pipe coverage layer and a response sequence for each subdivided segment is established. Then, the contribution weight sequence and phase difference sequence are obtained. Based on the contribution weight sequence and phase difference sequence, an exchange discriminant is constructed and compared with a threshold. The depth range and evidence indicators of the interlayer influence candidate segment are output. This combination of stratified contribution analysis and exchange evidence location makes the phase difference and mixed response differences primarily explained by the sequence of segment responses and interlayer exchange, rather than being converted into wellbore horizontal parameter anomalies. This mechanism reduces the probability of fitting the data but the layer correspondence not being valid. It also clearly presents the risks of interlayer crossflow, weak interlayer damage, and wellbore short circuit in the form of depth range.

[0017] Using candidate sections affected by interlayers as constraints, the vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity are jointly inverted. In-well observations are weighted and synthesized according to their contribution weight sequence for fitting. At the same time, the parameter degrees of freedom of non-critical sections are compressed through sharing rules to suppress overfitting. After the inversion is completed, the prediction error is recalculated under modified pumping conditions and compared with the error threshold. If the error exceeds the error threshold, a misfitting label is output. Thus, the output layer parameters correspond one-to-one with the depth layer, which can reduce the mismatch in the selection of curtain grouting depth, water-stopping layer, and recharge layer. Through threshold verification of the operating conditions extrapolation, the prediction reliability under different pumping conditions and different well network layouts is improved, and the engineering decision-making bias caused by extrapolation failure is reduced. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0019] Figure 1 This is a flowchart of a geostatistical inversion method for spatiotemporal heterogeneity of hydrogeological parameters according to the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] It should be noted that when a component is said to be "fixed to" another component, it can be directly attached to the other component or there may be an intervening component. When a component is said to be "connected to" another component, it can be directly connected to the other component or there may be an intervening component.

[0022] Please see Figure 1 As shown in the figure, this embodiment provides a geostatistical inversion method for spatiotemporal heterogeneity of hydrogeological parameters, including: S101. Obtain the water level time series and water quality time series of the monitoring well before pumping, during pumping, and during the recovery period. The water quality time series includes water quality time series such as conductivity and turbidity. Simultaneously obtain the pumping flow time series of the pumping well, as well as the top and bottom depths and interlayer thickness of the screen pipe covering layer of the monitoring well. Align the above data with a unified timestamp to form a hybrid observation dataset.

[0023] The method for obtaining water level time series and water quality time series such as conductivity and turbidity from monitoring wells before, during, and after pumping includes: setting up a water level detector at the wellhead and continuously recording the water level in the well at fixed sampling intervals, and continuously recording the water level segment before pumping to form a water level sequence before pumping; setting up conductivity and turbidity detectors on the sampling channel of the monitoring well, so that the conductivity and turbidity detectors continuously output conductivity and turbidity sequences from the start of pumping to the recovery period after pumping; recording the water level time series, conductivity sequence, and turbidity sequence output by the water level detectors using the same time stamp method to form a synchronous observation sequence of water level and water quality under the same pumping event; the reason is that retaining the time changes of water level response and water quality response simultaneously within the same pumping event provides input for establishing a time series correlation with the pumping flow rate time series.

[0024] The method for obtaining the pumping flow time series of a pumping well includes: installing a flow meter on the outlet pipeline of the pumping well to continuously output the instantaneous pumping flow and recording the flow changes corresponding to the start and end times of pumping; secondly, recording the output of the flow meter at the same sampling interval as the water quality time series such as turbidity as the pumping flow time series, and corresponding the pumping flow time series with the water level time series and the water quality time series using the same time identifier; describing the pumping process with quantifiable pumping disturbance intensity, so that the water level time series and the water quality time series such as conductivity and turbidity can be correlated with the pumping disturbance on the same time axis.

[0025] It should be noted that the pumping well serves as the point of application for pumping disturbance. The pumping flow rate time series is output according to the pumping process, so that the change in the pumping disturbance intensity over time can be quantified and traced. The pumping flow rate time series is used as an external excitation quantity on the time axis in the data organization. It is used to establish a temporal correlation between the water level time series, conductivity, turbidity and other water quality time series of the monitoring well and the same pumping event, so that the subsequent calculation of dynamic characteristics such as phase difference and hysteresis has a unified disturbance reference.

[0026] As observation points for the mixed response, monitoring wells, through the arrangement of screens spanning multiple aquifers or multiple subdivided segments, collect water level time series and water quality time series such as conductivity and turbidity, forming a mixed observation of the combined contribution of multiple layers after the propagation of pumping disturbances. At the same time, the monitoring wells provide the top and bottom depths of the screen-covered segments and the thickness of the interlayers, which are used to limit the vertical source range of the mixed observations. This allows subsequent processing to interpret the well observations as the result of a weighted synthesis of the screen-covered segments, rather than directly corresponding the observed values ​​to a single-layer simulated value, thereby reducing the risk of converting interlayer exchanges and the influence of weak interlayers into anomalies in horizontal parameters.

[0027] The combined role of pumping wells and monitoring wells is that pumping wells provide controllable and quantifiable disturbance inputs, while monitoring wells provide hybrid response outputs with vertical source constraints. After being aligned with a unified timestamp, they form a hybrid observation dataset, enabling subsequent steps to distinguish the differences in contributions from different layers and the phase differences caused by interlayer exchanges within the same pumping event framework, thus providing a basis for engineering risk identification and the selection of appropriate measures.

[0028] The method for obtaining the top and bottom depths and interlayer thicknesses of the screen pipe cover sections of monitoring wells includes: acquiring well completion data and well completion records of monitoring wells, and extracting the starting and ending depths of the screen pipes; overlaying the screen pipe depth range with the aquifer division information in the stratigraphic column or stratification record to obtain the top and bottom depths of each screen pipe cover section; identifying interlayers within the screen pipe cover area based on lithological records or well logging stratification results, and calculating the interlayer thickness based on the difference between the top and bottom depths of the interlayer; associating and storing the top and bottom depths of the screen pipe cover sections and the interlayer thicknesses with the pumping event identifiers used above, so that each pumping event corresponds to a set of screen pipe cover section information; and clarifying the vertical mixing source range corresponding to the water level time series and water quality time series such as conductivity and turbidity, providing spatial constraint information for the subsequent formation of a mixed observation dataset.

[0029] The method for forming a hybrid observation dataset includes: establishing a unified timestamp sequence using the time identifier of the pumping flow time series as a unified timestamp; mapping the water level time series and water quality time series such as conductivity and turbidity to the unified timestamp sequence, and filling in missing moments with linear interpolation of adjacent moments; calculating adjacent difference values ​​for abrupt changes in the well water level series and water quality series, and comparing the difference values ​​with a preset difference threshold. When the difference value exceeds the threshold, the data at that moment is replaced with the interpolation result of the adjacent moment; encapsulating the aligned well water level series, aligned conductivity series, aligned turbidity series, aligned pumping flow series, and the top and bottom depths and interlayer thickness of the screen cover section in the same recording unit to form a hybrid observation dataset; and solidifying the pumping disturbance intensity, well mixing response, and screen vertical cover information into the same data structure under a unified timestamp to provide common source input for subsequent stratified contribution decomposition.

[0030] Methods for calculating adjacent differences at abrupt changes in well water level and water quality sequences include: Based on a unified timestamp sequence, adjacent difference values ​​are calculated for the well water level sequence at two adjacent time points. The adjacent difference value is the difference between the water level measurements at adjacent times. Adjacent difference values ​​are also calculated for the conductivity and turbidity sequences at two adjacent time points to obtain water quality adjacent difference values. The adjacent difference values ​​for water level and water quality are compared with water level difference thresholds, and the adjacent difference values ​​for water quality are also compared with water quality difference thresholds. When any difference value exceeds the corresponding threshold, the corresponding time point is marked as a change point, and the timestamp and difference value of the change point are recorded. For the time marked as a change point, the measured value of the change point is replaced with the linear interpolation value of the adjacent time point, or the missing value is used in the mixed observation dataset. Loss markers are written for subsequent processing; after unified timestamp alignment, calculable adjacent change amplitudes are used to identify and suppress non-pumping mechanism abrupt changes caused by sensor instantaneous jitter, sampling disturbances, or instantaneous bubble entrainment, avoiding the abrupt change points being mistakenly regarded as evidence of interlayer exchange or interlayer influence when calculating phase difference increments, weighted abrupt change rates, and water level and water quality lags in subsequent calculations, thus ensuring the consistency of the time response characterization of the mixed observation dataset to pumping disturbances; for example, if the pumping flow rate changes smoothly but the turbidity suddenly spikes at a single timestamp, the difference value is replaced by interpolation after exceeding the threshold, which can avoid misjudging the spike as cross-layer transmission caused by particles released from interlayer channels.

[0031] S102. Based on the top and bottom depths of the screen tube cover layer, the vertical range corresponding to the mixed observation dataset is divided into subdivided layers, and the pumping flow time series is used as the excitation to establish a response sequence for each subdivided layer. The response sequences are matched with the water level time series and the water quality time series to obtain the contribution weight sequence and phase difference sequence of each subdivided layer.

[0032] Methods for dividing the vertical range corresponding to the mixed observation dataset into subdivided segments based on the top and bottom depths of the sieve tube cover layer include: Obtain the top and bottom depths of the sieve tube-covered strata, and determine the upper and lower boundary depths of the vertical range. Construct a strata boundary sequence based on the top and bottom depths of each stratum, and define basic strata between adjacent boundaries. Calculate the thickness of each basic stratum and compare it with a preset thickness threshold. Basic strata with thicknesses exceeding the thickness threshold are further subdivided into multiple sub-segments using equal depth steps. The equal depth step is obtained by rounding down the basic strata thickness and the required number of sub-segments. Generate a strata identifier for each sub-segment and record its top and bottom depths to form a set of sub-segments and associate it with the pumping event identifiers in the mixed observation dataset. Transform the vertical source range of the sieve tube coverage into a set of sub-segments with clear depth boundaries, so that subsequent contribution weights and phase differences can correspond to specific depth segments. For example, when the sieve tube covers three strata, three basic strata are obtained based on the top and bottom depths, and the basic strata with thicknesses exceeding the threshold are further subdivided into two segments using equal steps.

[0033] Methods for defining base segments between adjacent boundaries include: Obtain the segment boundary sequence and sort it from shallow to deep. The segment boundary sequence consists of the top and bottom depths of the sieve tube covered segment. Traverse the depth values ​​of two adjacent boundaries in the segment boundary sequence, take the shallower boundary depth as the top depth of the basic segment, and take the deeper boundary depth as the bottom depth of the basic segment, thus generating a basic segment record. Associate the basic segment record with the sieve tube covered segment identifier and write it into the basic segment table. The basic segment table contains at least the basic segment identifier, top depth, bottom depth, and thickness. If the top and bottom depths of the basic segment fall outside the starting and ending depths of the sieve tube, the basic segment record is discarded, and only the basic segment records that intersect with the vertical range of the sieve tube coverage are retained. Divide the vertical range of the sieve tube coverage into continuous and non-overlapping depth intervals with adjacent boundaries as the boundary, providing standardized segmentation units for further subdivision according to thickness thresholds. For example, if the boundary sequence is 10 meters, 18 meters, and 23 meters, then two basic segments are defined between 10 meters and 18 meters and between 18 meters and 23 meters, respectively.

[0034] Methods for calculating the thickness of each base layer include: The top and bottom depths of each basic segment in the basic segment table are obtained, and these depths are determined by the depths of adjacent boundaries. Next, the thickness of each basic segment is calculated by subtracting the top depth from the bottom depth, and the resulting thickness is associated with the corresponding basic segment identifier and written into the basic segment table. When the top and bottom depths use different depth benchmarks, they are converted to the same benchmark before calculating the basic segment thickness to ensure that the basic segment thickness is consistent with the top and bottom depths of the screen tube covering layer. The calculated basic segment thickness is compared with a preset thickness threshold, and the comparison result is written into the basic segment table for subsequent determination of whether the basic segment needs to be subdivided. This forms a basic segment thickness that can be directly used for threshold comparison, allowing the triggering conditions for subsequent subdivision of segments to be quantified and reproduced.

[0035] Methods for establishing response sequences for each sub-segment using pumping flow time series as excitation include: The pumping flow rate time series is acquired and uses the same timestamp sequence as the water level and water quality time series in the mixed observation dataset. For each sub-segment, the center depth and thickness of the segment are calculated based on the top and bottom depths. The center depth is mapped to a response delay, and the thickness is mapped to a response decay. The mapping uses a monotonic mapping rule and limits the range of the mapping results with preset delay and decay thresholds. A response kernel sequence for the sub-segment is constructed based on the response delay and decay, and the response kernel sequence is convolved with the pumping flow rate time series to obtain the response sequence of the sub-segment at the same timestamp. The above operation is repeated for the set of sub-segments to obtain a group of response sequences. The purpose is to convert the pumping flow rate time series into a computable response sequence for each sub-segment in the time domain, providing a common perturbation reference for subsequent matching with the water level and water quality time series. For example, under the same pumping flow rate input, sub-segments with shallower center depths are mapped to smaller delays, and their response sequences peak earlier.

[0036] Methods for calculating the center depth and thickness of a layer include: For each subdivided segment, the segment thickness is calculated by subtracting the segment top depth from the segment bottom depth, and the segment thickness is associated with the segment identifier of that subdivided segment and written into the segment parameter table. For each subdivided segment, the segment center depth is calculated by the arithmetic mean of the segment top depth and the segment bottom depth, and the segment center depth and segment thickness are written into the segment parameter table together. When the segment top depth and segment bottom depth use different depth references, the segment top depth and segment bottom depth are uniformly converted to the same depth reference before the calculation of the segment thickness and the segment center depth is performed to ensure consistency with the depth reference of the screen tube's starting and ending depths. The reason is that the calculable segment center depth and segment thickness characterize the spatial position and vertical scale of each subdivided segment, providing input for subsequently mapping the segment center depth to a response delay and the segment thickness to a response attenuation. For example, if a subdivided segment has a top depth of 20 meters and a bottom depth of 24 meters, then the segment thickness is 4 meters and the segment center depth is 22 meters.

[0037] Methods that employ monotonic mapping rules and limit the range of mapping results using preset delay and attenuation thresholds include: Obtain the center depth and thickness of the segment from the segment parameter table, and determine the center depth normalization interval based on the start and end depths of the screen tube, converting the segment center depth into a normalized value between zero and one. Generate an initial response delay value according to the monotonic mapping rule that a larger normalized center depth value corresponds to a larger response delay, and compare the initial response delay value with a preset delay threshold. If the response delay exceeds the threshold, set the response delay value to the preset threshold; otherwise, retain the initial response delay value as the response delay value corresponding to the segment center depth. Determine the thickness normalization interval based on the comparison result between the segment thickness and a thickness threshold, converting the segment thickness into a normalized value between zero and one. An initial value for response attenuation is generated according to a monotonic mapping rule that the larger the thickness normalization value, the greater the response attenuation. The initial value for response attenuation is compared with a preset attenuation threshold. If the response attenuation exceeds the threshold, the response attenuation is set to the preset attenuation threshold. If the response attenuation does not exceed the threshold, the initial value for response attenuation is retained, and the initial value for response attenuation is used as the response attenuation corresponding to the layer thickness. The response delay and the response attenuation are associated with the corresponding layer identifier and written into the layer parameter table as input for the subsequent construction of the response kernel sequence. For example, subdivided layers with a lower center depth receive a larger response delay, and subdivided layers with a larger thickness receive a larger response attenuation. When the mapping result exceeds the delay threshold or the attenuation threshold, it is restricted to the corresponding threshold range.

[0038] It should be noted that the initial values ​​generated by the monotonic mapping rule (i.e., the initial values ​​of response delay and response decay) have already transformed the propagation lag trend caused by center depth and the response decay trend caused by segment thickness into calculable parameters. When the initial values ​​fall within the threshold range (i.e., the delay threshold or decay threshold), it indicates that the parameter is within the physically reasonable range allowed by the threshold, and no truncation is necessary. Retaining the initial values ​​maintains the continuity and distinguishability of the mapping, ensuring that the response delay and response decay of different subdivided segments remain different, thereby facilitating the subsequent construction of the response kernel sequence. When matched with water level time series and water quality time series, it can reflect the differences in response time and response amplitude of each subdivision segment. If it is forcibly replaced with the threshold even when the threshold is not exceeded, the parameters of different segments will be compressed into the same value, which will weaken the distinguishability between subdivision segments and thus affect the decomposition results of subsequent contribution weight sequence and phase difference sequence. For example, when two segments have different center depths but do not exceed the delay threshold, retaining their initial values ​​can make their response peak times different. If they are both replaced with the same threshold, the two responses will converge in time and it will be difficult to distinguish the source of contribution.

[0039] Methods for constructing the response kernel sequence of this subdivision segment based on response delay and response decay include: Obtain the response delay and response decay associated with the subdivided segment identifier, and obtain the sampling interval of the unified timestamp sequence. Convert the response delay into a delay step, where the delay step is the integer part of the ratio of the response delay to the sampling interval. Fill the beginning of the response kernel sequence with zero values ​​according to the delay step, so that the response kernel sequence does not generate a response before the time corresponding to the delay step. Generate a decay segment after the delay step. The decay segment is generated according to a monotonically decreasing rule. The initial value of the decay segment is 1, and the subsequent values ​​gradually decrease according to the preset response decay, until the final value of the decay segment is not greater than the budgeted decay termination threshold or reaches the preset kernel length threshold. The process ends at the specified time. The response kernel sequence is normalized so that the sum of all terms in the response kernel sequence is one, and the response kernel sequence is associated with and stored with the subdivision segment identifier. The reason for this is that by fixing the response lag and response diffusion decay of the subdivision segment in the form of a discrete-time kernel, the subsequent effect on the pumping flow time series can be generated by a unified time-domain operation to produce candidate response trajectories for the subdivision segment. Moreover, the trajectory can be controlled in terms of time and amplitude by the response delay and response decay. For example, the response kernel sequence of a subdivision segment with a larger delay has a longer initial zero value, and the decay segment of a subdivision segment with a larger decay decreases faster.

[0040] The methods for setting the attenuation termination threshold and the kernel length threshold include: determining the duration window of disturbance impact that the response kernel needs to cover based on the sampling interval corresponding to the unified timestamp; setting candidate values ​​for the kernel length with reference to the change period of the pumping flow time series and the duration of the recovery segment of the water level time series and the water quality time series; and using the cumulative weight of the attenuation segment reaching a preset percentage threshold or the tail value of the attenuation segment being lower than the termination threshold as the termination condition, with the first of the two conditions being met triggering termination; and setting the two thresholds with the idea of ​​retaining the main response energy and suppressing tail noise.

[0041] Methods for performing convolution operations between the response kernel sequence and the pumping flow time series to obtain the response sequence of the subdivided segment under a unified timestamp include: Obtain the response kernel sequence associated with the subdivision segment identifier, and obtain the pumping flow time series associated with the same pumping event identifier and aligned with a unified timestamp sequence; wherein, the response kernel sequence consists of multiple kernel weights arranged in time step order, and the kernel weight is the value of the response kernel sequence at the corresponding position, used to characterize the contribution ratio of unit pumping disturbance at different backtracking time steps.

[0042] Confirm that the pumping flow time series and the response kernel series use the same sampling interval; when the sampling interval of the response kernel series is inconsistent with the sampling interval of the unified timestamp series, the response kernel series is resampled according to the sampling interval of the unified timestamp series, while keeping the structure of the response kernel series with zero delay segment and gradually decreasing attenuation segment unchanged.

[0043] For each target time in the unified timestamp sequence, a subdivided segment response value corresponding to that target time is generated. The generation method is as follows: the kernel weight of the response kernel sequence in the first position is multiplied by the pumping flow rate value of the target time; the kernel weight of the response kernel sequence in the second position is multiplied by the pumping flow rate value of the previous time step of the target time; and so on, following the order of the response kernel sequence, the pumping flow rate value of the corresponding time step back from the target time is multiplied by the kernel weight of the corresponding position; the multiplications obtained from each multiplication are accumulated and summed to obtain the subdivided segment response value of the target time, wherein the position number of the response kernel sequence corresponds one-to-one with the position number of the backtracking time step.

[0044] When the backward time step is earlier than the start time of the pumping flow time series, the pumping flow value of that time step is multiplied and accumulated as a zero value, so that it does not contribute before pumping begins. The multiplication and accumulation process is repeated for all target times of the unified timestamp sequence to obtain a convolutional output sequence covering the entire time period. The convolutional output sequence is then used as the subdivision segment response sequence under the unified timestamp sequence, and stored in association with the subdivision segment identifier and the pumping event identifier. An example of generating a response sequence is as follows: Given that the subdivision is a single segment, the sampling interval of the unified timestamp sequence is 10 seconds, the response delay of the subdivision is 30 seconds, the response decay of the subdivision is 0.2 (meaning that the decay segment decreases by 1 times the decay amount in each step), and the decay segment length is 5 steps (which can be determined by the decay termination condition or the kernel length control condition).

[0045] The delay steps are obtained from the response delay: Delay steps = response delay ÷ sampling interval = 30 ÷ 10 = 3 (rounded down).

[0046] Meaning: The first 3 sampling points of the response kernel sequence are 0, and no response is generated before the corresponding delay step.

[0047] Constructing an unnormalized response kernel sequence from the response decay: The attenuation factor r = 1 - response attenuation = 1 - 0.2 = 0.8. The attenuation segment (length 5 steps) is generated according to a monotonically decreasing rule. Step 1: 1; Step 2: 1 × 0.8 = 0.8; Step 3: 0.8 × 0.8 = 0.64; Step 4: 0.64 × 0.8 = 0.512; Step 5: 0.512 × 0.8 = 0.4096; Therefore, the unnormalized response kernel sequence h is: h=[0, 0, 0, 1, 0.8, 0.64, 0.512, 0.4096] Normalized response kernel sequence: Summation: S = 1 + 0.8 + 0.64 + 0.512 + 0.4096 = 3.3616; Normalized kernel sequence h′=h / S: h′=[0, 0, 0, 0.297477, 0.237982, 0.190386, 0.152308, 0.121847]; This means that after normalization, the sum of all terms is 1, making the magnitude of the convolution output comparable to the magnitude of the input. Subsequent magnitude threshold comparisons are only used as a constraint.

[0048] Given a pumping flow rate time series and performing convolution, the response sequence is obtained: Let the pumping flow rate time series x (the unit can be understood as cubic meters per hour) have 12 sampling points: x=[0, 0, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0] Corresponding time (seconds): 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110.

[0049] The convolution output y takes the same length as the causal output, i.e. , n is the nth sampling point, y[n] is the response value of the subdivision segment at the nth sampling point, x[nk] is the pumping flow input value at an earlier time (the nkth sampling point), h′[k] is the kernel weight of the response kernel sequence at the kth position, and Σ represents the summation of the effects of all historical inputs on the current time.

[0050] Calculation results (convolution completed): y=[0, 0, 0, 0, 0, 1.487387, 2.677297, 3.629224, 2.903379, 2.322703, 1.370776, 0.609234]; To clearly see the effect of delay and attenuation, key points are aligned and displayed: The pumping flow rate of 5 begins at t=20 seconds, but due to the delay step of 3, the response does not begin to rise until t=50 seconds.

[0051] t=50 seconds: y≈1.487 (corresponding to the first term 0.297477 of the input 5 after a delay when t=20 seconds enters the kernel: 5×0.297477≈1.487) t=60 seconds: y≈2.677 (combined with two inputs at t=20 and t=30: 5×0.237982+5×0.297477≈2.677) t=70 seconds: y≈3.629 (combining the three inputs t=20, t=30, and t=40: 5×0.190386+5×0.237982+5×0.297477≈3.629) The input is then zeroed out, but the tail of the core remains, so the response gradually decreases according to the decay phase.

[0052] The methods for matching each response sequence with water level time series and water quality time series to obtain the contribution weight sequence and phase difference sequence of each sub-segment include: Obtain the response sequence and water level time series for each sub-segment, and determine the peak time of the response sequence and the peak time of the water level time series for each sub-segment. A core interval is formed, starting from the minimum value of the peak time of the response sequence and the peak time of the water level time series for each sub-segment, and ending at the maximum value. The core interval is then extended by several time steps at both ends to obtain the main response time window. The number of extension steps is limited by a time window length threshold; steps exceeding the threshold are truncated to the threshold. This is to provide a calculation interval related to pumping disturbances for subsequent phase difference and weight correction. An example of the main response time window calculation is as follows: The sampling interval for the unified timestamp sequence is 10 seconds; the peak time of the response sequence for a certain subdivision segment is 70 seconds; and the peak time of the water level time sequence is 50 seconds.

[0053] The earliest peak time is 50 seconds, the later peak time is 70 seconds, and the core interval is [50 seconds, 70 seconds] (corresponding to sampling point numbers n=5 to n=7).

[0054] Let the two ends be extended by 2 time steps (i.e., 2 × 10 seconds = 20 seconds): starting point = 50 seconds - 20 seconds = 30 seconds, ending point = 70 seconds + 20 seconds = 90 seconds.

[0055] The main response window is [30 seconds, 90 seconds] (corresponding to n=3 to n=9, with a window length of 7 sampling points).

[0056] For each timestamp in the unified timestamp sequence, the response values ​​of each sub-segment at that timestamp are obtained and summed to obtain the segment response sum. The segment response sum is compared with a preset zero threshold. If the segment response sum is not greater than the zero threshold, the contribution weight of each segment at that timestamp is reset to zero. If the segment response sum is greater than the zero threshold, the response value of a certain sub-segment at that timestamp is divided by the segment response sum to obtain the contribution weight of that sub-segment at that timestamp. The contribution weight sequence of each sub-segment is output along the timestamps. The mixed response is directly allocated using the response proportion of each segment to obtain the contribution weight that changes over time. An example of contribution weight calculation is as follows: The sampling interval for the unified timestamp sequence is 10 seconds, and 5 timestamps are taken: t=0, 10, 20, 30, and 40 seconds; there are 3 subdivision segments: subdivision segment 1, subdivision segment 2, and subdivision segment 3; the response values ​​(sampled from the response sequence obtained by convolution) at each timestamp are shown in the table below:

[0057] Set the zero threshold to 0. Calculate by time stamp, t=0 seconds: the sum of the layer responses = 0 + 0 + 0 = 0. When it is not greater than the zero threshold, the contribution weights of the three sub-layers are all set to 0, and the contribution weights are: (0, 0, 0).

[0058] t=10 seconds: Segment response sum = 0.2 + 0.1 + 0.1 = 0.4, which is greater than the zero threshold, so the proportion is calculated: The weight of subdivided segment 1 is 0.2 / 0.4 = 0.50; The weight of subdivided segment 2 is 0.1 / 0.4 = 0.25; The weight of subdivided segment 3 is 0.1 / 0.4 = 0.25; Contribution weights: (0.50, 0.25, 0.25); t=20 seconds: Segment response sum = 0.6 + 0.3 + 0.1 = 1.0, which is greater than the zero threshold. The weight of subdivided segment 1 is 0.6 / 1.0 = 0.60; The weight of subdivided segment 2 is 0.3 / 1.0 = 0.30; The weight of subdivided segment 3 is 0.1 / 1.0 = 0.10; Contribution weights: (0.60, 0.30, 0.10); t=30 seconds: Segment response sum = 0.4 + 0.4 + 0.2 = 1.0, which is greater than the zero threshold: The weight of subdivided segment 1 is 0.4 / 1.0 = 0.40; The weight of subdivided segment 2 is 0.4 / 1.0 = 0.40; The weight of subdivided segment 3 is 0.2 / 1.0 = 0.20; Contribution weights: (0.40, 0.40, 0.20); t=40 seconds: When the sum of the responses of the three sub-segments is 0+0+0=0 and is not greater than the zero threshold, the contribution weights of the three sub-segments are all set to 0, and the contribution weights are (0, 0, 0). The output contribution weight sequence is arranged by timestamp: Subdivision 1: [0, 0.50, 0.60, 0.40, 0], Subdivision 2: [0, 0.25, 0.30, 0.40, 0], Subdivision 3: [0, 0.25, 0.10, 0.20, 0].

[0059] For each sub-segment, the peak moment of the response sequence of that sub-segment is calculated within its main response time window, and the peak moment of the water level time sequence within the same main response time window is also calculated. The difference between the two peak moments is converted into a step difference, which is the rounded ratio of the time difference to the sampling interval. The above calculation of the peak moment difference is repeated along the time axis by a sliding time window, and the phase difference sequence of that sub-segment is output and stored in association with the sub-segment identifier. The reason for this is to quantify the chronological relationship between the segment response and the water level response using the peak moment difference, forming a directly comparable phase difference sequence. An example of phase difference sequence calculation is as follows: The sampling interval for the unified timestamp sequence is 10 seconds, and the timestamps are taken as t=0, 10, 20, 30, 40, 50, 60, 70, 80, and 90 seconds. The response sequence of a certain subdivision segment (corresponding to the above timestamps) is: R=[0, 0, 1, 3, 5, 4, 2, 1, 0, 0]. The water level time series is: W=[0, 1, 2, 4, 6, 5, 3, 1, 0, 0]. The sliding window length is 5 sampling points (i.e., 50 seconds), and the sliding step size is 1 sampling point (i.e., 10 seconds).

[0060] Time window 1: Coverage t = 0–40 seconds (sampling points 0 to 4); R subsequence: [0, 0, 1, 3, 5], peak value is 5, peak time = 40 seconds; W subsequence: [0, 1, 2, 4, 6], peak value is 6, peak time = 40 seconds; Time difference = 40 - 40 = 0 seconds; Step difference = (0 ÷ 10) rounded up = 0 steps; Phase difference value: 0 steps.

[0061] Time window 2: Coverage t = 10-50 seconds (sampling points 1 to 5); R subsequence: [0, 1, 3, 5, 4], peak value is 5, peak time = 40 seconds; W subsequence: [1, 2, 4, 6, 5], peak value is 6, peak time = 40 seconds; Time difference = 0 seconds, step difference = 0 steps, phase difference value: 0 steps.

[0062] Time window 3: Coverage t = 20-60 seconds (sampling points 2 to 6); R subsequence: [1, 3, 5, 4, 2], peak value = 5, peak time = 40 seconds; W subsequence: [2, 4, 6, 5, 3], peak value = 6, peak time = 40 seconds; Time difference = 0 seconds, step difference = 0 steps, phase difference value: 0 steps.

[0063] Time window 4: Coverage t = 30–70 seconds (sampling points 3 to 7); R subsequence: [3, 5, 4, 2, 1], peak value = 5, peak time = 40 seconds; W subsequence: [4, 6, 5, 3, 1], peak value = 6, peak time = 40 seconds; Time difference = 0 seconds, step difference = 0 steps, phase difference value: 0 steps.

[0064] Time window 5: Coverage t = 40-80 seconds (sampling points 4 to 8); R subsequence: [5, 4, 2, 1, 0], peak value = 5, peak time = 40 seconds; W subsequence: [6, 5, 3, 1, 0], peak value = 6, peak time = 40 seconds; Time difference = 0 seconds, step difference = 0 steps, phase difference value: 0 steps.

[0065] This yields the phase difference sequence for this subdivided segment (output according to the start timestamp of each time window): Phase difference sequence = [0, 0, 0, 0, 0] (unit: steps) Examples of controls with lag are as follows: Maintaining a sampling interval of 10 seconds, the overall peak value of the response sequence in this subdivided segment lags the peak water level by 20 seconds: R2 = [0, 0, 0, 1, 3, 5, 4, 2, 1, 0] (peak 5 at 50 seconds); W remains: [0, 1, 2, 4, 6, 5, 3, 1, 0, 0] (peak 6 at 40 seconds); within a time window covering t = 20-60 seconds: R2 subsequence (20-60 seconds) = [0, 1, 3, 5, 4], peak time = 50 seconds; W subsequence (20-60 seconds) = [2, 4, 6, 5, 3], peak time = 40 seconds; Time difference = 50 - 40 = +10 seconds The step difference = (10 ÷ 10) rounded up = +1 step, the phase difference is +1 step, which means that the response of this section lags behind the water level by 1 sampling interval.

[0066] Obtain the conductivity and turbidity sequences, calculate adjacent difference sequences, and compare the difference values ​​with a preset water quality difference threshold. Timestamps where the difference value exceeds the threshold are marked as water quality change moments. For each water quality change moment, obtain the phase difference sequence of each sub-segment, calculate the corresponding response arrival time, and calculate the absolute value of the time difference between the response arrival time and the water quality change moment. Compare the absolute value of the time difference with the time difference threshold; segments whose absolute value does not exceed the threshold are designated as matching segments. Adjust the contribution weight sequence at the water quality change moment by increasing the weight of matching segments according to a preset weight adjustment step size, and proportionally decreasing it for non-matching segments, ensuring the sum of contribution weights for all segments remains one at that moment. This is because, while maintaining the response proportion sharing framework, the contribution weight is tilted towards time-aligned segments at the moment of water quality change, thus keeping the weight sequence consistent with the water quality change in time. An example of weight adjustment for matching segments is shown below: The unified timestamp sampling interval is 10 seconds, with t=40, 50, 60, 70, and 80 seconds; the two water quality sequences are as follows: Conductivity sequence C(t): [100, 102, 140, 141, 142]; Turbidity sequence T(t): [5.0, 5.1, 5.0, 5.0, 5.1]; The water quality differential thresholds are set as follows: conductivity differential threshold = 30, turbidity differential threshold = 1; at t = 60 seconds, the phase difference between the two sub-layer segments (unit: steps) is: Subdivision segment 1: Phase difference = 0 steps; Subdivision segment 2: Phase difference = +2 steps.

[0067] The time difference threshold is set to 15 seconds; the weight correction step size is set to 0.10; at this moment (t=60 seconds), the contribution weights before correction are: segment 1=0.60, segment 2=0.40 (the sum of the two is 1).

[0068] Calculate the adjacent differences and mark the time of water quality change: Calculate the adjacent difference ΔC for the conductivity sequence: t=50: ΔC=102-100=2, which does not exceed the conductivity difference threshold; t=60: ΔC=140-102=38, which exceeds the conductivity difference threshold of 30. Mark t=60 as the moment of water quality change.

[0069] t=70: ΔC=141-140=1, which does not exceed the conductivity difference threshold; t=80: ΔC=142-141=1, which does not exceed the conductivity difference threshold.

[0070] Calculate the adjacent difference ΔT for the turbidity sequence: t=50: ΔT=5.1-5.0=0.1, which does not exceed the turbidity difference threshold; t=60: ΔT=5.0-5.1=-0.1, which does not exceed the turbidity difference threshold. ... Therefore, the water quality change time in this example is 60 seconds.

[0071] Calculate the arrival time of the response from the phase difference and compare the absolute values ​​of the time differences: Response arrival time = Water quality change time + Phase difference × Sampling interval); Subdivision 1: Phase difference = 0 steps; Response arrival time = 60 + 0 × 10 = 60 seconds; The absolute value of the time difference = |60 - 60| = 0 seconds; Compare with the time difference threshold of 15 seconds: if 0 does not exceed the time difference threshold of 15 seconds, then subdivision segment 1 is the matching segment; Subdivision 2: Phase difference = +2 steps Response arrival time = 60 + 2 × 10 = 80 seconds; The absolute value of the time difference = |80-60| = 20 seconds; Comparison with the time difference threshold of 15 seconds: 20 exceeds the time difference threshold of 151, so subdivision segment 2 is a non-matching segment.

[0072] Adjust the contribution weights upward or downward at the moment of water quality change, while keeping the sum equal to one: Weights before correction at t=60 seconds: w1=0.60, w2=0.40.

[0073] Adjust the weighted segment 1 upwards by a step size of 0.10: After the adjustment, the weight w1′=0.60+0.10=0.70.

[0074] For the unmatched segment 2, the weights are proportionally reduced to keep the sum of the weights equal to one: Since there is only one unmatched segment remaining, we directly set w2′=1-w1′=0.30. The corrected weights are w1′=0.70 and w2′=0.30, and the sum of the two is still 1.

[0075] S103. Calculate the exchange discriminant based on the contribution weight sequence and the phase difference sequence. The exchange discriminant includes the phase difference increment, weight mutation rate, water level and water quality lag. Compare the exchange discriminant with the corresponding preset threshold. When any discriminant exceeds the corresponding preset threshold, output the depth range of the interlayer influence candidate segment.

[0076] The sampling interval and full-time length of the unified timestamp sequence are obtained, and a time window index is generated according to the sliding step size. Each time window is determined by the start timestamp and the end timestamp. The time window length and sliding step size are compared with the time window length threshold and the step size threshold, respectively, and the values ​​that do not exceed the threshold are taken. The set of subdivided segments and their top depth and bottom depth are obtained, and the contribution weight sequence and phase difference sequence associated with each subdivided segment identifier are obtained. At the same time, the water level time series and water quality time series associated with the same pumping event identifier are obtained. The time window index is associated with the timestamp, subdivided segment, contribution weight and phase difference to form an alignment table. The reason for this is that each subsequent exchange discrimination quantity can be calculated and traced within the determined time window and determined subdivided segment.

[0077] Methods for calculating phase difference increment and weight mutation rate based on contribution weight sequence and phase difference sequence include: The subdivided segments are sorted from shallow to deep according to the center depth of the segment, and adjacent subdivided segment pairs are defined. Within each time window, the phase difference sequence values ​​of adjacent subdivided segment pairs are obtained, the phase difference increment is calculated, and the phase difference increment is the phase difference value of the lower subdivided segment minus the phase difference value of the upper subdivided segment. The phase difference increment is then associated with the boundary depth of the two segments. Within the same time window, the contribution weight sequence of each sub-segment is obtained. The absolute value of the difference in contribution weight between adjacent timestamps within the time window is calculated, and the absolute value of the difference is divided by the sampling interval to obtain the weight mutation rate. This forms a weight mutation rate sequence, which is then correlated with the depth range of the sub-segment. The reason for this is that the temporal mutation of the response difference between adjacent segments and the contribution distribution of segments is quantified into two components of the exchange discriminant, which facilitates subsequent threshold comparison to locate the influence of the interlayer. An example of calculating the phase difference increment and the weight mutation rate is as follows: The sampling interval Δt = 10 seconds. The subdivision segments are sorted from shallowest to deepest according to their center depth: segment A (upper part), segment B (lower part), and the boundary depth between them is 18 meters. Therefore, the adjacent subdivision segment pairs are (A, B). The phase difference increment will be associated with the 18-meter boundary depth. A time window W2 is taken: t = 20 seconds to t = 80 seconds (including timestamps 20, 30, 40, 50, 60, 70, 80). Within this time window, the phase difference sequence values ​​(unit: steps; 1 step = 10 seconds) and the contribution weight sequence values ​​are given: Phase difference sequence (P): t=20: P A =1, P B =3; t=40: P A =2, P_B=4; t=60: P A =2, P B =5; t=80: P A =3, P B =6.

[0078] Contribution weight sequence (W): t=20: W A =0.70, W B =0.30; t=40: W A =0.60, W B =0.40; t=60: W A =0.40, W B =0.60; t=80: W A =0.30, W B =0.70.

[0079] Calculate the phase difference increment within time window W2: Phase difference increment definition: ΔP AB(t) =P B(t)- P A(t) Time-by-time calculation: t=20: ΔP AB =3-1=2 steps (corresponding to 20 seconds); t=40: ΔP AB =4-2=2 steps (corresponding to 20 seconds); t=60: ΔP AB =5-2=3 steps (corresponding to 30 seconds); t=80: ΔP AB =6-3=3 steps (corresponding to 30 seconds); therefore, the phase difference increment sequence is obtained within the time window W2: ΔP AB =[2, 2, 3, 3], in steps, and store the sequence in association with a boundary depth of 18 meters.

[0080] For example, if the phase difference of the lower segment is 2 to 3 steps greater than that of the upper segment, it indicates that the lower segment is more lagging behind the upper segment, and the lag difference is concentrated near the 18-meter boundary.

[0081] Calculate the weighted mutation rate within time window W2: Adjacent timestamps use 20, 40, 60 to 80; For segment A: 20 to 40: |0.60 - 0.70| ÷ 20 = 0.10 ÷ 20 = 0.005 / second; 40 to 60: |0.40 - 0.60| ÷ 20 = 0.20 ÷ 20 = 0.010 / second; 60 to 80: |0.30-0.40|÷20=0.10÷20=0.005 / second.

[0082] The weighted mutation rate sequence of segment A within W2 is obtained: R WA =[0.005, 0.010, 0.005] / second, and associated with the depth range of 10–18 meters in section A.

[0083] For segment B: 20 to 40: |0.40 - 0.30| ÷ 20 = 0.10 ÷ 20 = 0.005 / second; 40 to 60: |0.60 - 0.40| ÷ 20 = 0.20 ÷ 20 = 0.010 / second; 60 to 80: |0.70 - 0.60| ÷ 20 = 0.10 ÷ 20 = 0.005 / second; The weighted mutation rate sequence of segment B within W2 is obtained: R WB =[0.005, 0.010, 0.005] / second, and associated with the depth range of 18–23 meters in section B.

[0084] For example, in the 40 to 60 second range, the weight mutation rate of both layers reaches 0.010 / second, indicating that the contribution distribution of the layers has been significantly adjusted within this time window. Subsequently, it can be compared with the weight mutation rate threshold to trigger the interlayer influence candidate.

[0085] Methods for calculating water level and water quality lags based on water level time series and water quality time series include: Within each time window, adjacent difference sequences are calculated for the water level time series, and adjacent difference sequences are calculated for the conductivity and turbidity sequences in the water quality time series. Cross-correlation values ​​are calculated between the water level and conductivity difference sequences, and the time step count corresponding to the maximum cross-correlation value is taken as the lag for water level and conductivity. Similarly, cross-correlation values ​​are calculated between the water level and turbidity difference sequences, and the time step count corresponding to the maximum cross-correlation value is taken as the lag for water level and turbidity. These lags are associated with the time window identifier and stored to form the water level and water quality lags. The purpose is to quantify the sequential relationship between water level and water quality changes using the alignment method of difference changes within the same time window, allowing the lags to be directly used as the third term in the exchange discrimination metric for threshold comparison. An example of calculating the water level and water quality lags is shown below. The sampling interval Δt = 10 seconds; the time window identifier W2 covers t = 20, 30, 40, 50, 60; the upper limit threshold for the number of time shift steps is 1 step, so the number of cross-correlation time shift steps m takes {-1, 0, +1}.

[0086] Water level time series H=[100, 98, 96, 95, 95]; The conductivity sequence is C = [500, 500, 505, 510, 512]. The turbidity sequence T = [10, 10, 10, 12, 14].

[0087] Calculate the adjacent difference sequence, according to the definition of adjacent difference ΔX[i]=X[i]-X[i-1], within the time window W2: ΔH=[-2,-2,-1,0]; ΔC=[0,5,5,2]; ΔT=[0,0,2,2].

[0088] Calculate the cross-correlation sequence, and for each pair of difference sequences, calculate the cross-correlation in the range m∈{-1, 0, +1}: R HC(m =∑ΔH[i]×ΔC[i+m] (outbound terms are not included in the summation); R HT(m) =∑ΔH[i]×ΔT[i+m] (outbound terms are not included in the summation).

[0089] The calculation yielded: R HC(0) =-15, R HC(+1 )=-22,R HC(-1) =-5; R HT(0) =-2, R HT(+1) =-6, R HT(-1) =0.

[0090] The number of time steps corresponding to the peak cross-correlation is taken as the hysteresis. For R... HC(m) With R HT(m) Calculate the cross-correlation magnitude |R (m) | Select the m with the largest amplitude as the number of time steps; this m is the step representation of the lag: For conductivity: |R HC(+1) |=22 is the maximum, water level and conductivity hysteresis = +1 step = 10 seconds; for turbidity: |R HT(+1) |=6 is the maximum, and the lag between water level and turbidity is +1 step = 10 seconds.

[0091] The hysteresis sequence is formed by storing the time window identifier together. The result of W2 is recorded as follows: W2, water level and conductivity hysteresis = +1 step, water level and turbidity hysteresis = +1 step.

[0092] Other time windows are processed in the same way and written into records to obtain the water level and water quality lag sequence.

[0093] When any discriminant exceeds the corresponding preset threshold, the method for outputting the depth range of the candidate segment affected by the mezzanine includes: The phase difference increment, weighted mutation rate, and water level and water quality lag are obtained for each time window. The phase difference increment is compared with the phase difference increment threshold, the weighted mutation rate is compared with the weighted mutation rate threshold, and the water level and water quality lag is compared with the lag threshold. When any discriminant exceeds the corresponding preset threshold, the type of discriminant triggered, the discriminant value, the time window identifier, and the associated sub-segment identifier or the boundary depth of adjacent segments are recorded to form an exchange overthreshold set. The reason is that the triggering evidence of the interlayer effect is generated by a unified threshold comparison rule to avoid the use of relative judgment.

[0094] For each record in the exchange threshold set, if the trigger type is phase difference increment (i.e., exceeding the phase difference increment threshold), then the adjacent subdivision segment pairs associated with the record are obtained, and the top and bottom depths of the two segments are merged to obtain the depth range of the candidate segment, centered on the boundary depth of the two segments; if the trigger type is weighted mutation rate (i.e., exceeding the weighted mutation rate threshold), then the subdivision segment identifier associated with the record is obtained, and the top and bottom depths of the subdivision segment are output as the depth range of the candidate segment; if the trigger type is water level and water quality lag... If the hysteresis threshold is exceeded, the identifier of the segment with the highest contributing weight sequence within that time window is obtained, and the top and bottom depths of that segment are output as the depth range of the candidate segment. The depth range and the corresponding discriminant value are output together and stored in association with the pumping event identifier. The reason is to convert the time-domain evidence of exceeding the threshold into a depth range in the vertical space, and use the discriminant value as traceable evidence to provide input for subsequent vertical inversion constrained by the candidate segment. An example of outputting the depth range and discriminant value of the interlayer affecting the candidate segment is as follows: Subdivisions: Section A has a top depth of 10 meters and a bottom depth of 18 meters; Section B has a top depth of 18 meters and a bottom depth of 23 meters; Section C has a top depth of 23 meters and a bottom depth of 30 meters; the boundary depths between adjacent sections are 18 meters and 23 meters; window index: time window W2 is 20–80 seconds; time window W3 is 40–100 seconds.

[0095] Exchange discrimination parameters within each time window: W2: Phase difference increment ΔP AB =2 steps, ΔP BC =1 step; weighted mutation rate R WA =0.004 / second, R WB =0.012 / second, R WC =0.003 / second; hysteresis L HC =+1 step, L HT =+1 step.

[0096] W3: Phase difference increment ΔP AB =1 step, ΔP BC =3 steps; weighted mutation rate R WA =0.006 / second, R WB =0.007 / second, R WC =0.010 / second; hysteresis L HC =+3 steps, L HT =+2 steps.

[0097] The preset phase difference increment threshold is 2 steps; the preset weight mutation rate threshold is 0.010 / second; the preset hysteresis threshold is 2 steps; the judgment rule is that if the discriminant is greater than the corresponding threshold, it is considered an over-threshold.

[0098] Threshold alignment and formation of exchanged superthreshold sets: For time window W2: Phase difference increment: ΔP AB =2 steps without exceeding the threshold; ΔP BC =1 step without exceeding the threshold.

[0099] Weighted mutation rate: R WB =0.012 / second exceeding the threshold, record one exchange threshold exceeding record: R1 = {Pumping event identifier E1, time window identifier W2, trigger type weight mutation rate, discriminant value 0.012 / second, subdivision segment identifier B}.

[0100] Lag: L HC =+1 step does not exceed the threshold; L HT =+1 step without exceeding the threshold.

[0101] Regarding time window W3: Phase difference increment: ΔP BC =3 steps to exceed the threshold, record one exchange-based threshold exceedance record: R2 = {Pumping event identifier E1, time window identifier W3, trigger type phase difference increment, discrimination value taking 3 steps, adjacent subdivision segment pair (B, C), boundary depth 23 meters}.

[0102] Weighted mutation rate: R WC =0.010 / second does not exceed the threshold.

[0103] Lag: L HC =+3 steps exceeding the threshold; to determine the segment with the maximum weight, take the average contribution weight of each segment within the time window W3 as W. A Mean 0.30, W B Mean 0.55, W C With a mean of 0.15 and the maximum weight segment identified as B, record one instance of exchange exceeding the threshold: R3 = {pumping event identifier E1, time window identifier W3, trigger type water level and water quality lag, discrimination value L} HC =+3 steps, maximum weight segment identifier B}.

[0104] Therefore, the exchanged overthreshold set in this example contains three records R1, R2, and R3.

[0105] Output the depth range of candidate segments affected by the mezzanine based on the trigger type: R1, whose weighted mutation rate exceeds the threshold, outputs the depth range of layer B as 18–23 meters, and outputs a discriminant value of 0.012 / second, forming an evidence record: E1 / W2: Candidate segment 18–23 meters; Evidence is a weighted mutation rate of 0.012 / second. It should be noted that the E1 / W2 positioning index is used to bind the evidence to a specific pumping event and a specific time window. The same applies to E1 / W2 and E1 / W3 below.

[0106] R2, whose phase difference increment exceeds the threshold, outputs the depth range of the adjacent layer pair (B, C) after merging. The merged range is from 18 meters at the top of B to 30 meters at the bottom of C, i.e., 18–30 meters. It also outputs the discriminant value in 3 steps and the boundary depth of 23 meters, forming an evidence record. E1 / W3: Candidate section 18–30 meters; evidence is 3-step phase difference increment; boundary depth 23 meters.

[0107] R3, whose water level and water quality lag exceed the threshold, outputs the maximum weighted layer B at a depth range of 18–23 meters.

[0108] E1 / W3: Candidate section 18–23 meters.

[0109] S104. Establish a vertical inversion model with the depth range of the candidate section affected by the interlayer as a constraint. Combine the water level time series and water quality time series according to the contribution weight sequence to form a weighted observation, and form a weighted observation sequence set to participate in the fitting. Invert the vertical permeability, interlayer conductivity and wellbore short-circuit conductivity. Under the condition of changing the pumping flow rate time series, recalculate the prediction error and compare the prediction error with the preset error threshold to identify the misaligned fitting.

[0110] Methods for establishing a vertical inversion model using the depth range of candidate segments affected by interlayers as constraints include: Obtain the depth range of interlayer influence candidate segments associated with the same pumping event identifier, and obtain the top and bottom depths of the subdivided segment sets; mark subdivided segments that intersect with the depth range as candidate segments, and mark subdivided segments that do not intersect as non-candidate segments; establish a layered parameter set for the vertical inversion model, so that vertical permeability is independently valued within candidate segments and according to a preset sharing rule within non-candidate segments; limit interlayer conductivity to independently value only adjacent segments that cross the boundary depth of candidate segments, and other adjacent segments are valued according to a sharing rule; set parameter boundary threshold constraints for wellbore short-circuit conductivity, and associate this parameter with the top and bottom depths of the screen pipe covered segment. The purpose is to use the candidate segment depth range to constrain the parameter degree of freedom allocation, so that the inversion focus is on the relevant segments and boundaries of the candidate segments, and avoid introducing unnecessary independent parameters in non-candidate segments. A specific example is as follows: Subdivided layer sets: S1 (top depth 10 m, bottom depth 14 m), S2 (top depth 14 m, bottom depth 18 m), S3 (top depth 18 m, bottom depth 23 m), S4 (top depth 23 m, bottom depth 30 m). The boundary depths of adjacent layers are 14 m, 18 m, and 23 m; the top and bottom depths of the sieve tube covered layers are: sieve tube starting depth 10 m, sieve tube ending depth 30 m; the depth range of candidate zones affected by interlayers is 18–23 m; the intersection thickness threshold is δ, taken as δ = 0.1 m, and the judgment rule is: if the intersection thickness between a subdivided layer and a candidate zone is greater than δ, it is judged as having an intersection; otherwise, it is judged as having no intersection.

[0111] Candidate and non-candidate segments are marked by depth intersection, and the intersection thickness of each subdivision segment and the candidate segment (18–23 meters) is calculated: S1 intersects with the candidate segment with a thickness of 0 meters, not greater than δ, and is marked as a non-candidate segment; S2 intersects with the candidate segment with a thickness of 0 meters, not greater than δ, and is marked as a non-candidate segment; S3 intersects with the candidate segment with a thickness of 5 meters, greater than δ, and is marked as a candidate segment; S4 intersects with the candidate segment with a thickness of 0 meters, not greater than δ, and is marked as a non-candidate segment.

[0112] Therefore, the candidate segment set is {S3}, and the non-candidate segment set is {S1, S2, S4}.

[0113] Establish a hierarchical parameter set for the vertical inversion model and determine the independent or shared values ​​of vertical permeability: Candidate segment S3: Set vertical permeability as a segment-independent parameter Kz. S3 Non-candidate segment layers are merged into shared groups according to preset sharing rules; the sharing rules use two threshold comparisons: the thickness difference threshold is Δh. th The center depth difference threshold is Δz th Adjacent layers must simultaneously satisfy the condition that the thickness difference does not exceed Δh. th And the difference in center depth does not exceed Δz th They are then merged into the same shared group. For ease of demonstration, let Δh be taken. th =1 meter, Δz th =5 meters; Calculate adjacent layer pairs in non-candidate sections: S1 and S2: thicknesses are 4 meters and 4 meters respectively, with a thickness difference of 0 meters not exceeding Δh. th The center depths are 12 meters and 16 meters respectively, and the difference in center depth is 4 meters, which does not exceed Δz_th. Therefore, S1 and S2 are merged into a shared group G1, sharing the parameter Kz. G1 And assign values ​​to S1 and S2 simultaneously; S2 and S4 are not adjacent segment pairs. S3 exists in between and is not included in the adjacent merging judgment. S4 forms a separate shared group G2, sharing the parameter Kz. G2And assign the value to S4.

[0114] Therefore, the parameter to be inverted for vertical permeability is Kz. S3 Kz G1 Kz G2 Kz G1 Simultaneously representing S1 and S2.

[0115] The independent value of the inter-layer conductivity is determined at certain locations, while a shared value is used for the remaining locations. The adjacent layer pairs are (S1, S2), (S2, S3), and (S3, S4), corresponding to boundary depths of 14 meters, 18 meters, and 23 meters, respectively. The boundary depths of the candidate segment 18–23 meters are 18 meters and 23 meters. Therefore, an independent parameter G for the inter-layer conductivity is set for the adjacent layer pairs (S2, S3) that cross the 18-meter boundary. 18 The subscript 18 signifies the boundary of a segment with a depth of 18 meters. This means the boundary depth between adjacent segments (S2, S3) is 18 meters. Therefore, the inter-layer conductivity at this boundary is denoted as G. 18 For adjacent layer pairs (S3, S4) that cross a 23-meter boundary, an independent inter-layer conductivity number G is set. 23 For adjacent layer pairs (S1, S2) that do not cross the boundary of candidate segments, set a shared parameter G. out .

[0116] Therefore, the parameter to be inverted for interlayer conduction capability is G. 18 G 23 G out .

[0117] To set a parameter boundary threshold constraint for the wellbore short-circuit conductivity, and to correlate it with the top and bottom depths of the screen pipe cover layer, the wellbore short-circuit conductivity is set as parameter G. well And apply parameter boundary threshold constraints to it, the constraint form being: G well With lower bound threshold G welllow and upper bound threshold G wellhigh Comparison, making G well Not less than G welllow and not greater than G wellhigh At the same time, G well The effective depth range is limited to between 10 meters at the beginning of the screen and 30 meters at the end of the screen, so that the model only takes into account the conductivity within the range of 10–30 meters when calculating the short-circuit related flux in the wellbore.

[0118] In this example, the hierarchical parameter set of the vertical inversion model is: Vertical permeability: Kz S3 (Candidate segment independent), Kz G1 (Shared by S1 and S2), Kz G2(S4 shared); Interlayer connectivity: G 18 (18-meter boundary independent), G 23 (23-meter boundary independent), G out (14-meter shared boundary); Shaft short-circuit conductivity: G well (subject to G) welllow With G wellhigh The constraint limits the effective depth to 10–30 meters within the screen tube.

[0119] Therefore, the independently valued parameters only appear in candidate segment S3 and its two adjacent boundaries at depths of 18 meters and 23 meters. The other non-candidate segments and boundaries adopt shared values, thereby limiting the independent parameters to the relevant segments and boundaries of the candidate segments and reducing the number of independent parameters in non-candidate segments.

[0120] It should be noted that the methods for setting the values ​​of the preset sharing rules include: Obtain the water level and water quality time series corresponding to each pumping event from the historical pumping event identifier set, and obtain the associated subdivided segment sets, segment top depth and segment bottom depth, contribution weight sequence and phase difference sequence; calculate the thickness difference and center depth difference of adjacent segment pairs based on the subdivided segment sets, and compare the thickness difference with the thickness difference threshold and the center depth difference with the depth difference threshold one by one, and count the number of adjacent segment pairs that are merged into the same shared group in each pumping event to form a merging frequency table; iteratively adjust the thickness difference threshold and depth difference threshold according to the merging frequency table to ensure that the prediction error after merging in each pumping event does not exceed the error threshold. The proportion of shared groups reaches the proportion threshold, and the number of shared groups is not less than the number of groups threshold. The adjusted thickness difference threshold and depth difference threshold are solidified as preset sharing rule values ​​and associated with the top and bottom depth range of the screen tube coverage layer and the depth range of the candidate section and written into the rule table. In a new pumping event, the adjacent layer pairs of non-candidate sections are compared and merged according to the thickness difference threshold and depth difference threshold in the rule table to obtain the shared group identifier and drive the assignment of shared parameters of non-candidate sections. For example, if under a certain threshold, the prediction error does not exceed the error threshold in eight out of ten pumping events after merging and the number of shared groups meets the number of groups threshold, then the threshold of that group is used as the preset sharing rule value.

[0121] The method of combining water level time series and water quality time series according to contribution weights to form a weighted observation set, and using this weighted observation set to participate in fitting, and inverting vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity includes: Obtain the contribution weight sequence associated with each sub-segment identifier, and obtain the water level time series and water quality time series associated with the same pumping event identifier and aligned with a unified timestamp. For each timestamp, weight the water level observation, conductivity observation, and turbidity observation values ​​according to the contribution weight of each sub-segment to form a weighted observation sequence set. The weighted observation sequence set includes a weighted water level observation sequence, a weighted conductivity observation sequence, and a weighted turbidity observation sequence, and the weighting process keeps the sum of the contribution weights of each segment equal to one. Input the pumping flow time series into the vertical inversion model, and output a prediction sequence group aligned with the unified timestamp sequence. The prediction sequence group includes at least a water level prediction sequence, a conductivity prediction sequence, and a turbidity prediction sequence. Calculate the water level prediction sequence... The difference sequences between the column and weighted water level observation sequences are calculated, as are the difference sequences between the conductivity prediction sequence and the weighted conductivity observation sequence, and the difference sequences between the turbidity prediction sequence and the weighted turbidity observation sequence. The sum of the squares of each difference sequence is then combined to form a fitting error index. Under parameter boundary threshold constraints, the vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity are iteratively updated, causing the fitting error index to decrease successively until it no longer decreases or reaches the iteration number threshold. This is because the contribution weights are used to convert the mixed well observations into weighted observations, allowing both water level and water quality observations to jointly constrain the vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity. This distinguishes the differences in the contributions of parameters inside and outside the candidate section to the error term. A specific example is as follows: The unified timestamps are t0, t1, and t2; the water level time series is H=[10, 9, 8]; the conductivity series is C=[100, 110, 130]; the turbidity series is T=[5, 6, 8]; the contribution weight sequences of the two subdivided segments A and B are wA=[0.7, 0.5, 0.3] and wB=[0.3, 0.5, 0.7], and each timestamp satisfies wA+wB=1. The pumping flow time series Q=[0, 5, 5] is used as the excitation input for the vertical inversion model.

[0122] Weighted observations are synthesized based on the contribution weight sequence to form a weighted observation sequence set. For each timestamp, the observation value is multiplied by the segment weight to obtain the weighted observation component sequence. The entire component sequence set is then defined as the weighted observation sequence. Weighted water level observation sequence: Weighted water level observation sequence H of section A A =[7.0, 4.5, 2.4], weighted water level observation sequence H of section B B =[3.0, 4.5, 5.6]; Weighted conductivity observation sequence: Weighted conductivity observation sequence C for layer A A =[70, 55, 39], Weighted conductivity observation sequence C of segment B B =[30, 55, 91]; Weighted turbidity observation sequence: Weighted turbidity observation sequence T of segment A A =[3.5, 3.0, 2.4], weighted turbidity observation sequence T of segment B B =[1.5, 3.0, 5.6].

[0123] The above weighting process is performed separately for water level observations, conductivity observations, and turbidity observations, with each timestamp having a weight sum of one.

[0124] The vertical inversion model is input with the pumping flow rate time series, and outputs a set of predicted sequences aligned with a unified timestamp sequence. Water level prediction sequence: Water level prediction sequence for section A =[6.6, 4.7, 2.6], Water level prediction sequence for section B =[3.2, 4.3, 5.4]; Conductivity prediction sequence: Conductivity prediction sequence for segment A =[72, 54, 41], Predicted conductivity sequence for segment B =[28, 56, 89]; Turbidity prediction sequence: Turbidity prediction sequence for segment A =[3.6, 2.9, 2.5], Turbidity prediction sequence for segment B =[1.4, 3.1, 5.5].

[0125] The sum of squares of each difference sequence is then combined to form a fitting error index: Water level difference in section A The square of the water level in section A and SE HA =0.24; water level difference in section B The square of the water level in section B and SE HB =0.12; Subtotal of water level SE H =SE HA +SE HB =0.36.

[0126] Conductivity difference of segment A The sum of squares of the conductivity of segment A, SE CA =9; Conductivity difference of segment B The sum of squares of the conductivity of segment B, SE CB =9; Subtotal of conductivity SE C =SE CA +SE CB =18.

[0127] Turbidity difference in section A: The sum of squares of turbidity SE in section A TA =0.03; Turbidity difference in section B The sum of squares of turbidity SE in section B TB =0.03; Turbidity subtotal SE T =SE TA +SE TB =0.06.

[0128] The fitting error index is the sum of three factors: SE total =SE H +SE C +SE T =18.42.

[0129] Under parameter boundary threshold constraints, the vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity are iteratively updated, and the fitting error index is recalculated. Parameter update step sizes are set and limited by parameter step size thresholds. Update candidates for vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity are updated by the current value ± the step size. Each update candidate is compared with its corresponding lower and upper bound thresholds; if it exceeds the bounds, the boundary value is used. The updated parameters are input into the vertical inversion model to recalculate the prediction sequence group, and the SE is recalculated. total If SE is recalculated total Smaller than the previous SE total If the update is successful, the update is retained and the process proceeds to the next round; otherwise, the original parameters are retained and the iteration count is incremented. When SE total The process stops when the iteration count no longer decreases or reaches the iteration threshold, and outputs the final vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity. The method for recalculating the prediction error under varying pumping flow rate time series conditions and comparing the prediction error with a preset error threshold to identify misfitting layers includes: The modified pumping flow time series conditions are constructed as follows: The pumping flow time series corresponding to another pumping event identifier is obtained, or the amplitude transformation and time shift are performed on the pumping flow time series under the same pumping event identifier to generate the modified pumping flow time series. While maintaining the vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity obtained from the inversion, the modified pumping flow time series is input into the vertical inversion model to obtain the recalculated water level prediction series and the recalculated water quality prediction series. The recalculated water level prediction series and the corresponding water level time series under the corresponding conditions are then calculated respectively. The difference sequence of the column is used to calculate the difference sequence between the recalculated predicted water quality sequence and the water quality time series under the corresponding conditions, and the sum of squares of the difference sequence is combined into the prediction error. The prediction error is compared with a preset error threshold. When the prediction error exceeds the preset error threshold, a misfitting indicator is output, and the depth range of the candidate segment that triggers the indicator and the corresponding stratification parameter values ​​are associated with the output. The purpose is to use changes in pumping flow conditions to form excitation inputs different from the original fit, so that parameter combinations that fit only a single pumping condition are identified as misfitting when the conditions change due to the prediction error exceeding the threshold. A specific example is as follows: Candidate section depth range: 18–23 meters, associated with pumping event identifier E1; The values ​​of the inverted vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity parameters are recorded in P1 (remaining unchanged during recalculation): Vertical permeability Kz S3 =0.8、Kz G1 =1.2、Kz G2 =0.9; Interlayer conductivity G 18 =0.5, G 23 =0.4, G out =0.2; Wellbore short-circuit conductivity G well =0.3; Uniform timestamps: t0, t1, t2, t3; Preset error threshold: E th =12.

[0130] K represents the permeability parameter, z represents the vertical direction; S3 represents the subdivision segment identifier S3, that is, the subdivision segment corresponding to a specific depth range; Kz S3 This indicates the vertical penetration capability of the subdivided segment; G1 and G2 represent the shared group identifiers, obtained by merging several non-candidate segment subdivisions using sharing rules; Kz G1 This indicates the vertical permeability shared by all sub-segments within the shared group G1; G represents the inter-layer conductivity, and 18 and 23 represent the boundaries between adjacent segments with boundary depths of 18 meters and 23 meters, respectively; G 18 G indicates the inter-layer conductivity at the 18-meter junction. 23 Similarly; out represents the non-candidate segment boundary sharing term, used for sharing inter-layer connectivity between adjacent layer pairs that do not cross the candidate segment boundary depth; G wellThis is a parameter for the short-circuit conductivity of the wellbore. "well" indicates that the conductivity corresponds to the short-circuit channel in the wellbore, rather than the interlayer interface.

[0131] Give the pumping flow time series and observation series of the original pumping event identifier E1, and calculate the original condition prediction error.

[0132] Pumping flow rate time series Q for pumping event identifier E1 E1 =[0, 5, 5, 0]; the observation sequence associated with E1 is the water level time series H. E1 =[10, 9, 8, 8], conductivity sequence C E1 =[100, 110, 130, 128], Turbidity sequence T E1 =[5, 6, 8, 7]; Keep parameter value record P1 unchanged, and set Q... E1 Input the vertical inversion model to obtain the water level prediction sequence. Conductivity prediction sequence Turbidity prediction sequence Specifically =[10.0, 9.1, 8.1, 7.9], =[100, 112, 128, 127], =[5.0, 6.1, 7.9, 7.2]; Calculate the sum of squared differences: SE H for The sum of the squares of each point yields 0.03; SE C for The sum of the squares of each point is 9; SE T for The sum of the squares obtained point by point is 0.06; the sum of the three values ​​gives the prediction error E of the pumping event identifier E1. E1 =9.09, and compared with the preset error threshold E th =12 comparison, E E1 It does not exceed the threshold.

[0133] Construct the modified pumping flow time series conditions and obtain the corresponding observation series: Using another pumping event identifier E2 as the changing condition: pumping flow rate time series Q E2 =[0, 0, 7, 7]; the observation sequence associated with E2 is H E2 =[10, 10, 9, 8]、C E2 =[100, 100, 108, 120]、T E2 =[5, 5, 5.5, 6.5], and with Q E2 Align on t0 to t3.

[0134] Keeping P1 unchanged, Q E2Input the vertical inversion model and obtain =[10.0, 9.8, 8.4, 7.6], =[100, 103, 120, 140], =[5.0, 5.4, 7.2, 9.0].

[0135] Calculate SE H2 for The sum of squares at each point is 0.56; calculate SE. C2 for The sum of squares obtained point by point is 553; calculate SE. T2 for The sum of the squares obtained point by point is 9.30; adding the three together gives the prediction error E of the pumping event identifier E2. E2 =562.86. E E2 With E th =12 comparison, E E2 Exceeding the threshold triggers misaligned fitting.

[0136] The output record is fixed as follows: {pumping event identifier = E2, misaligned fitting identifier = 1, prediction error = E} E2 Preset error threshold = E th The candidate segment depth range is 18–23 meters, and the stratification parameter value record is P1}, where P1 includes Kz. S3 Kz G1 Kz G2 G 18 G 23 G out With G well The value of .

[0137] It should be noted that the prediction error originates from the sum of squares of the differences between the recalculated water level prediction sequence and the corresponding water level time series, as well as the sum of squares of the differences between the recalculated water quality prediction sequence and the corresponding water quality time series. The recalculated water level prediction sequence and the recalculated water quality prediction sequence are calculated by the vertical inversion model after inputting the changed pumping flow rate time series. The corresponding water level time series and the corresponding water quality time series are derived from a mixed observation dataset and aligned with a unified timestamp. Weighted observations are used to weight and synthesize the water level time series and the water quality time series according to the contribution weight sequence. Therefore, the prediction error is data-dependent on the pumping flow rate time series, the water level time series, the water quality time series, the contribution weight sequence, the depth range of the candidate section affected by the interlayer, and the vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity obtained from the inversion. For example, if the contribution weight sequence is consistently high in the corresponding layer of the candidate section, and the inversion parameters result in an excessively high value for the wellbore short-circuit conductivity, then the sum of squares of the water quality prediction difference will accumulate and increase under the condition of changing the pumping flow rate time series, leading to the prediction error exceeding the threshold.

[0138] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters, characterized in that, include: S101. Obtain the water level time series and water quality time series of the monitoring well before pumping, during pumping, and during the recovery period. Simultaneously obtain the pumping flow time series of the pumping well, as well as the top and bottom depths and interlayer thickness of the screen pipe covering layer of the monitoring well. Align the above data with a unified timestamp to form a hybrid observation dataset. S102. Based on the top and bottom depths of the screen tube cover layer, the vertical range corresponding to the mixed observation dataset is divided into subdivided layers, and the pumping flow time series is used as the excitation to establish a response sequence for each subdivided layer; each response sequence is matched with the water level time series and the water quality time series to obtain the contribution weight sequence and phase difference sequence of each subdivided layer. S103. Calculate the exchange discriminant based on the contribution weight sequence and the phase difference sequence; The exchange discrimination value is compared with the corresponding preset threshold. When any discrimination value exceeds the corresponding preset threshold, the depth range of the interlayer influence candidate segment is output. S104. Establish a vertical inversion model with the depth range of the candidate section affected by the interlayer as a constraint. Combine the water level time series and water quality time series according to the contribution weight sequence to form a weighted observation, form a weighted observation sequence set to participate in the fitting, and invert the vertical permeability, interlayer conductivity and well short-circuit conductivity. Under the condition of changing the pumping flow rate time series, recalculate the prediction error and compare the prediction error with the preset error threshold to identify the misaligned fitting.

2. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 1, characterized in that, Methods for dividing the vertical range corresponding to the mixed observation dataset into subdivided segments based on the top and bottom depths of the sieve tube cover segment include: Obtain the top and bottom depths of the sieve tube covering layer to determine the upper and lower boundary depths of the vertical range, and construct the layer boundary sequence based on the top and bottom depths of each layer. Obtain the segment boundary sequence and sort it by depth. Traverse the two adjacent boundary depth values ​​in the segment boundary sequence, take the shallower boundary depth as the top depth of the basic segment and take the deeper boundary depth as the bottom depth of the basic segment to generate the basic segment table, and remove the basic segments that do not intersect with the start and end depths of the sieve tube. The top and bottom depths of the basic layer segments are obtained from the basic layer segment table, and the thickness of the basic layer segment is calculated by subtracting the top depth from the bottom depth. The thickness of the basic layer segment is then compared with a preset thickness threshold. Basic layers whose thickness exceeds the thickness threshold are subdivided into multiple sub-segments according to equal depth steps.

3. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 2, characterized in that, Methods for establishing response sequences for each sub-segment using pumping flow time series as excitation include: Obtain the pumping flow rate time series and use a unified timestamp sequence; For each subdivided segment, the center depth and thickness of the segment are calculated based on the top and bottom depths of the segment. The center depth of the segment is mapped to a response delay according to a monotonic mapping rule, and the mapping result is compared with a delay threshold to limit the mapping result. The thickness of the segment is mapped to a response decay according to a monotonic mapping rule, and the mapping result is compared with a decay threshold to limit the mapping result. The response kernel sequence of the subdivided segment is constructed based on the response delay and response decay. The response kernel sequence is convolved with the pumping flow time series to obtain the response sequence of the subdivided segments.

4. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 3, characterized in that, Methods for obtaining subdivided response sequences by performing convolution operations on the response kernel sequence and the pumping flow time series include: Obtain the response kernel sequence associated with the subdivision segment identifier, and obtain the pumping flow time series associated with the same pumping event identifier and aligned with a unified timestamp sequence; The sampling interval of the response kernel sequence is compared with the sampling interval of the unified timestamp sequence. If the sampling intervals are inconsistent, the response kernel sequence is resampled according to the sampling interval of the unified timestamp sequence, while keeping the delay segment at zero and the decay segment gradually decreasing. For each target moment in the unified timestamp sequence, the pumping flow rate value at the corresponding time step is traced back according to the position of the response kernel sequence and multiplied with the corresponding kernel weight. The results are then accumulated to obtain the subdivided segment response value. When the traced back time step is earlier than the start time of the pumping flow rate time series, the multiplication and accumulation are performed with a zero value. The response kernel sequence consists of multiple kernel weights arranged in time step order. The kernel weight is the value of the response kernel sequence at the corresponding position and is used to characterize the contribution ratio of a unit pumping disturbance at different traced back time steps. The response sequence of each subdivided segment is obtained by repeatedly multiplying and accumulating the responses at all target times.

5. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 3, characterized in that, Methods for obtaining the contribution weight sequence and phase difference sequence of each sub-segment include: Obtain the response sequence and water level time series of each sub-segment, determine the peak time of the response sequence and the peak time of the water level time series of each sub-segment, and determine the core interval by taking the minimum value of the peak time as the starting point and the maximum value of the peak time as the ending point. Based on the core interval, a main response window is generated at both ends by an extended time step. Within the main response window, the response values ​​of each subdivided segment are summed according to a unified timestamp to obtain the segment response sum. The contribution weight is determined based on the comparison between the sum of the segment response and the zero threshold. When the sum of the segment response is not greater than the zero threshold, the contribution weight of each sub-segment at the corresponding timestamp is reset to zero. When the sum of the segment response is greater than the zero threshold, the sub-segment response value at the corresponding timestamp is divided by the sum of the segment response to obtain the sub-segment contribution weight at that timestamp. The contribution weight sequence of each sub-segment is output along the unified timestamp. Based on the main response time window, the peak time of the response sequence of the subdivided segment and the peak time of the water level time sequence are determined respectively. The difference between the two peak times and the sampling interval are converted into the step difference. The phase difference sequence is obtained by repeatedly converting along the time axis by sliding time window. The water quality time series includes conductivity and turbidity sequences. The conductivity and turbidity sequences are obtained and adjacent difference sequences are calculated separately. The difference values ​​are compared with the water quality difference threshold. Timestamps exceeding the water quality difference threshold are marked as water quality change times. At each water quality change time, the phase difference sequence values ​​are read to calculate the absolute value of the time difference between the response arrival time and the water quality change time. The absolute value of the time difference is compared with the time difference threshold to determine the matching segment. At the water quality change time, the matching segment is adjusted upward according to the preset weight, and the non-matching segment is adjusted downward proportionally, while keeping the sum of the contribution weights of each sub-segment equal to one.

6. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 1, characterized in that, Methods for calculating the exchange discriminant based on contribution weight sequences and phase difference sequences include: The exchange discrimination parameters include phase difference increment, weighted abrupt rate, and water level and water quality lag. Obtain the sampling interval and full time period length of the unified timestamp sequence, generate a time window index according to the sliding step size, and determine each time window by the start timestamp and the end timestamp; Obtain the subdivided layer set and the top and bottom depths of each layer. Also obtain the contribution weight sequence and phase difference sequence associated with each subdivided layer identifier. Obtain the water level time series and water quality time series associated with the same pumping event identifier. Establish an alignment table by associating the time window index with the timestamp, subdivided layer, contribution weight, and phase difference. The subdivided layer set is composed of subdivided layers. Based on the alignment table, adjacent subdivided segments are sorted by the center depth of the segment to form pairs. The phase difference increment is calculated in each time window and correlated with the boundary depth of the two segments. The absolute value of the difference between adjacent timestamps of the contribution weight sequence is calculated and divided by the sampling interval to obtain the weight mutation rate sequence, which is then correlated with the depth range of the subdivided segment. Based on the alignment table, the adjacent difference sequences of the water level time series and the adjacent difference sequences of the water quality time series are calculated in each time window. The cross-correlation value sequence is calculated and the time step number corresponding to the maximum cross-correlation value is taken as the lag between water level and conductivity and the lag between water level and turbidity. These are associated with the time window identifier and stored to form the lag between water level and water quality.

7. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 6, characterized in that, When any discriminant exceeds the corresponding preset threshold, the method for outputting the depth range of the candidate segment affected by the mezzanine includes: The phase difference increment, weighted mutation rate, and water level and water quality lag are obtained in each time window. The phase difference increment is compared with the phase difference increment threshold, the weighted mutation rate is compared with the weighted mutation rate threshold, and the water level and water quality lag is compared with the lag threshold. When any discriminant exceeds the corresponding threshold, the discriminant type, discriminant value, time window identifier, and associated subdivision segment identifier or adjacent segment boundary depth are recorded to form an exchange overthreshold set. For each record in the exchange threshold set, when the discriminant type is phase difference increment, the top depth and bottom depth of the two segments are merged based on the boundary depth of adjacent subdivision segments to obtain the candidate segment depth range. When the discriminant type is weighted mutation rate, the top depth and bottom depth of the subdivision segment are output based on the subdivision segment identifier to obtain the candidate segment depth range. When the discriminant type is water level and water quality lag, the top depth and bottom depth of the segment with the largest weight in the time window contribution weight sequence are output based on the water level and bottom depth of the segment to obtain the candidate segment depth range.

8. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 1, characterized in that, Methods for establishing a vertical inversion model using the depth range of candidate segments affected by interlayers as constraints include: Obtain the depth range of interlayer impact candidate segments associated with the same pumping event identifier, and obtain the top depth and bottom depth of the subdivided segment set. Based on the intersection relationship between the depth range and the subdivided segments, mark the candidate segment segments and non-candidate segment segments. A layered parameter set for the vertical inversion model is established based on candidate and non-candidate segments. The layered parameter set includes vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity. Vertical permeability is taken independently by segment within the candidate segment and by a preset sharing rule within the non-candidate segment. Based on the boundary depth of the candidate segment, the adjacent layer pairs that cross the boundary depth of the candidate segment are determined. The inter-layer connectivity is set independently for the adjacent layer pairs and set according to the preset sharing rules for the other adjacent layer pairs. Set parameter boundary threshold constraints for wellbore short-circuit conductivity and correlate wellbore short-circuit conductivity with the top and bottom depths of the screen pipe cover section.

9. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 8, characterized in that, The method of combining water level time series and water quality time series according to contribution weights to form a weighted observation set, and using this weighted observation set to participate in fitting, and inverting vertical permeability, interlayer conductivity, and wellbore short-circuit conductivity includes: Obtain the contribution weight sequence associated with each sub-segment identifier and obtain the water level time series and water quality time series associated with the same pumping event identifier and aligned with a unified timestamp; For each timestamp, the water level observation, conductivity observation and turbidity observation are weighted and synthesized according to the contribution weight of each sub-segment to form a weighted observation sequence set, and the weighting process keeps the sum of the contribution weights of each segment equal to one. The vertical inversion model is input with the pumping flow time series and outputs a set of prediction sequences aligned with a unified timestamp sequence, including at least the water level prediction sequence, conductivity prediction sequence and turbidity prediction sequence. The difference sequences between the water level prediction sequence and the weighted water level observation sequence, the difference sequences between the conductivity prediction sequence and the weighted conductivity observation sequence, and the difference sequences between the turbidity prediction sequence and the weighted turbidity observation sequence are calculated separately, and the sum of the squares of each difference sequence is combined as the fitting error index. Under the constraint of parameter boundary threshold, the vertical permeability, interlayer conductivity and wellbore short-circuit conductivity are iteratively updated, and the fitting error index is recalculated until the fitting error index no longer decreases or reaches the iteration number threshold. The vertical permeability, interlayer conductivity and wellbore short-circuit conductivity are then output.

10. The method for geostatistical inversion of spatiotemporal heterogeneity of hydrogeological parameters according to claim 9, characterized in that, The method for recalculating the prediction error under varying pumping flow rate time series conditions and comparing the prediction error with a preset error threshold to identify misfitting layers includes: The modified pumping flow time series is constructed by obtaining the pumping flow time series corresponding to another pumping event identifier, or by performing amplitude transformation and time shift on the pumping flow time series under the same pumping event identifier to generate the modified pumping flow time series. While keeping the vertical permeability, interlayer conductivity and wellbore short-circuit conductivity unchanged, the modified pumping flow rate time series is input into the vertical inversion model to output the recalculated water level prediction series and the recalculated water quality prediction series. Calculate the difference sequence between the recalculated water level prediction sequence and the water level time series under the corresponding conditions, and calculate the difference sequence between the recalculated water quality prediction sequence and the water quality time series under the corresponding conditions. Combine the sum of squares of the difference sequences to form the prediction error. The prediction error is compared with a preset error threshold. When the prediction error exceeds the preset error threshold, a misfitting identifier is output. The candidate segment depth range that triggers the misfitting identifier and the corresponding layer parameter values ​​are output and associated and stored.

Images