Intelligent control method and system for mine grouting equipment
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2026-06-12
- Publication Date
- 2026-08-11
AI Technical Summary
[0007]为解决现有矿井注浆控制方法在地质工况变化时无法有效感知数据分布漂移及其成因、且难以将浆液渗流规律等物理机理约束融入模型更新的缺陷,导致控制适应性不足的问题,本发明提供了一种矿井注浆设备智能控制方法及系统
通过并行计算监测数据的高阶统计量变化率与模型特征SHAP归因值时序梯度,构建对工况变化具有高灵敏度的多维漂移特征向量,并以马氏距离精确判定分布偏移,当偏移超限时,由漂移模式分类器识别具体地质成因并定位待更新功能模块,再从预置模型库中实例化符合当前流固关联规律的物理约束方程,将其融入条件变分自编码器生成机理驱动的虚拟注浆样本。后续采用基于欧氏距离的负指数加权策略,为不同关联度的历史样本差异化分配训练权重,结合漂移程度动态调整学习率,仅对受影响模块执行局部参数更新。该方法有效解决了传统方案对早期渐变型漂移检测不灵敏、更新过程忽略物理规律、小样本下模型适配滞后及全量重训引发模型历史有效知识遗忘的问题,显著提升了复杂矿井地质条件下注浆控制的准确性与响应及时性。
Smart Images

Figure CN122386725B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of intelligent control technology in mining engineering. More specifically, this invention relates to an intelligent control method and system for mine grouting equipment. Background Technology
[0002] During mine grouting operations, as the mining face advances, geological conditions such as the degree of rock fissure development, permeability, and stress distribution can change continuously or suddenly. This places extremely high demands on the control strategies of grouting equipment. If the control model cannot detect and respond to these changes in a timely manner, key process parameters such as grouting pressure and grout diffusion range are prone to deviating from safe ranges, potentially leading to grouting failure or even water inrush accidents.
[0003] To address the aforementioned issues, several adaptive control schemes for grouting under varying working conditions have been proposed in the prior art. For example, Chinese patent document CN112127918B discloses an automatic grouting control method and device. This scheme is applicable to mine roof grouting for water inrush prevention. It integrates multiple execution units and various types of sensors to collect real-time parameters of the grouting process, and adopts a segmented threshold control strategy to control the grouting volume and pressure according to the two different stages of initial grouting and pressure injection. After the grout density reaches a set range in the initial grouting stage, it transitions to the pressure injection stage. By monitoring the change trend of the orifice pressure, the grouting effect and the timing of stage transition are determined, ultimately achieving segmented adjustment of the grouting process parameters.
[0004] However, the existing technical solutions mentioned above still have the following shortcomings when facing the grouting scenario in mines where geological conditions are constantly changing: First, the parameter adjustment strategy of the existing methods is essentially a direct response to the deviation of sensor measurements from preset thresholds. It does not detect whether the working condition characteristics have shifted from the data distribution level. When the geological conditions in the mining area change gradually or suddenly, the statistical distribution characteristics of the monitoring data and the degree of dependence of the control model on each input feature will change. The existing methods can only passively wait for the measured value to exceed the threshold before triggering the control action. They cannot capture the implicit working condition drift trend in time, nor can they identify the cause type of drift event, resulting in a delayed and untargeted control response.
[0005] Second, the adaptive capability of existing methods is entirely based on preset stage thresholds and switching rules. When the working conditions change beyond the preset range, it is impossible to introduce physical constraints such as the seepage equation of grout in the surrounding rock fissures and mass conservation. Therefore, it is difficult to generate virtual grouting samples that conform to the characteristics of the new working conditions to assist in model updates. When integrating historical data to adjust the model, the correlation differences between different historical data and the characteristics of the current working conditions are not considered, making it difficult to ensure the engineering rationality of control parameters under the new working conditions.
[0006] In summary, there is an urgent need for an intelligent control scheme for mine grouting equipment that can effectively sense data distribution drift and identify its causes under varying working conditions, while incorporating the inherent physical constraints of the grouting process, such as slurry seepage patterns, into the control model update process, in order to solve the problem of insufficient adaptability of existing technologies under changing geological conditions. Summary of the Invention
[0007] To address the shortcomings of existing mine grouting control methods, such as the inability to effectively perceive data distribution drift and its causes when geological conditions change, and the difficulty in incorporating physical mechanisms such as grout seepage patterns into model updates, resulting in insufficient control adaptability, this invention provides an intelligent control method and system for mine grouting equipment.
[0008] This invention provides an intelligent control method for mine grouting equipment, comprising: S1, acquiring a multidimensional online monitoring data stream of the mine grouting process, and outputting grouting equipment control commands from an initial control model; generating a multidimensional drift feature vector based on the rate of change of higher-order statistics of the multidimensional online monitoring data stream within a sliding time window and the time-series gradient of the SHAP attribution value of the input features of the initial control model; S2, when the first distance between the multidimensional drift feature vector and the historical stable state feature distribution exceeds a preset threshold, identifying the cause of the drift event and the functional modules to be updated in the initial control model through a drift pattern classifier, and updating the initial control model according to the cause. S3. Obtain the appropriate constraint equation from the model library; S4. Integrate the constraint equation as a loss term into the data generation model to generate virtual grouting process data, and output the corresponding virtual control command from the initial control model; S5. Construct an enhanced training set by fusing the virtual grouting process data, the corresponding virtual control command, and the pre-acquired historical monitoring data; S6. Assign non-uniform training weights based on the second distance between the current multidimensional drift feature vector and the multidimensional drift feature vector corresponding to the historical monitoring data, adjust the learning rate in combination with the first distance to locally update the functional module to be updated, and output the grouting equipment control command from the updated control model.
[0009] A multidimensional drift feature vector is constructed by fusing the rate of change of higher-order statistics from grouting monitoring data with the temporal gradient of the SHAP attribution value of model input features. Mahalanobis distance is used to determine geological condition drift, accurately identifying the causes of drift and the modules to be updated in the model. The fluid-solid correlation constraint equation adapted to the new working conditions is incorporated as a loss term into the conditional variational autoencoder to generate virtual grouting data that conforms to physical laws to supplement the new working condition samples. Non-uniform training weights are assigned to historical samples based on Euclidean distance, and the learning rate is dynamically adjusted by combining Mahalanobis distance and only locally updated for the target module. This can shorten the model update cycle, reduce the model's forgetting of historical working conditions, and improve the stability and adaptability of grouting control under complex mine geology.
[0010] Preferably, the method based on the rate of change of higher-order statistics of the multidimensional online monitoring data stream within a sliding time window includes: setting the time step and window length of the sliding time window, and acquiring the multidimensional online monitoring data stream of each stage in real time according to the time step; calculating the skewness and kurtosis values of each dimension of the multidimensional online monitoring data stream within the current sliding time window; using a backward difference algorithm to calculate the difference between the skewness and kurtosis values of each corresponding dimension between the current sliding time window and the immediately preceding sliding time window; dividing each difference by the time step to obtain the rate of change of skewness and kurtosis of each dimension of the monitoring data corresponding to each sliding time window, which is used as the rate of change of higher-order statistics.
[0011] By setting the time step and window length of a sliding time window, the multidimensional online monitoring data stream of mine grouting is segmented in real time. The skewness and kurtosis values of each dimension of data within the current window are calculated. Then, the difference between the corresponding statistics of adjacent windows is calculated using a backward difference algorithm and divided by the time step to obtain the skewness and kurtosis change rates of each dimension of monitoring data as higher-order statistical change rates. This method can characterize the drift of the centroid of grouting condition data distribution and the abrupt change rate of the distribution tail extrema, and can more sensitively capture early gradual geological condition changes, providing reliable statistical characteristics for the accurate determination of subsequent conceptual drift.
[0012] Preferably, the temporal gradient of the SHAP attribution value of the input feature of the initial control model is obtained through the following steps: An analytical algorithm is used to interpret the inference phase of the initial control model, and an independent baseline SHAP attribution value is calculated for each feature dimension of the input feature within the sliding time window; using the baseline SHAP attribution value corresponding to each feature dimension within the current sliding time window, the baseline SHAP attribution value of the corresponding feature dimension in the immediately preceding sliding time window is subtracted to obtain the first-order discrete difference; the first-order discrete differences of all feature dimensions are divided by the time step of the sliding time window and combined into a vector, which serves as the temporal gradient of the SHAP attribution value of the input feature.
[0013] Preferably, the step of identifying the causes of drift events and the functional modules to be updated in the initial control model through the drift pattern classifier includes: inputting the multidimensional drift feature vector into the drift pattern classifier, calculating the predicted probability vector components of each category label generated by the output layer; extracting the category index node where the component with the largest value is located; reading the dictionary hash table with pre-configured static association information, querying through the category index node mapping, and obtaining the specific flow field structure variation state features as the cause of drift; and obtaining the corresponding independent network layer parameter names by matching the field columns corresponding to the cause, so as to determine the functional modules to be updated.
[0014] This method combines data-driven pattern classification with causal mapping guided by prior knowledge. After detecting a shift in working condition distribution, it can automatically trace the anomaly back to the specific geomechanical root cause and accurately determine the functional components in the control model that need to be adjusted, thus avoiding the computational cost of a comprehensive parameter update.
[0015] Preferably, the pre-set model library is a fluid-structure interaction model library. The step of identifying the causes of drift events and the functional modules to be updated in the initial control model through a drift pattern classifier, and obtaining suitable constraint equations from the pre-set model library according to the causes, includes: extracting partial differential flow field parameter equations with matching names and containing unassigned initial coefficient terms from the fluid-structure interaction model library according to the identified causes; extracting the slurry viscosity sequence and instantaneous flow rate sequence from the current multidimensional online monitoring data stream, and obtaining the three-dimensional spatial coordinates corresponding to each grouting hole location and pipeline sensor; substituting the temporal characteristics of the slurry viscosity sequence into the fluid dynamic viscosity coefficient term of the partial differential flow field parameter equation, and substituting the instantaneous flow rate sequence into the dynamic flow boundary condition term of the corresponding spatial coordinate point, completing parameter assignment, and generating the constraint equations through discretization processing and spatial topology mapping.
[0016] Preferably, the generation of the constraint equations through discretization and spatial topology mapping includes: calling a mathematical discretization algorithm component to transform the continuous partial differential flow field parameter equations with initial coefficients into a set of linear algebraic equilibrium equations between grid nodes according to the divided spatial grid; constructing a dimension-reduced mapping observation operator matrix based on the three-dimensional spatial coordinates of each grouting hole location and pipeline sensor to form a set of constraint equations matching the current geological conditions as the constraint equations.
[0017] By integrating the mechanistic model with the field sensing topology, the generated constraint equations can not only express the inherent laws of slurry seepage and surrounding rock deformation, but also be strictly aligned with the actual observable data in the spatial dimension, providing a constraint benchmark that conforms to physical reality for the subsequent generation of virtual data.
[0018] Preferably, the data generation model is a conditional variational autoencoder. Integrating the constraint equations as loss terms into the data generation model to generate virtual grouting process data includes: inputting the multidimensional online monitoring data stream and the multidimensional drift feature vector into the conditional variational autoencoder to generate a virtual full-space grid state tensor and the mean and variance of the latent variable distribution; subsequently, substituting the virtual full-space grid state tensor into the constraint equations to obtain the mechanism penalty function term, and extracting data from it using the observation operator matrix to calculate the native reconstruction loss; finally, adding the mechanism penalty function term, the native reconstruction loss, and the KL divergence term derived based on the mean and variance to obtain the total loss; updating the network parameters of the conditional variational autoencoder based on the total loss, and generating the virtual grouting process data based on random noise and the multidimensional drift feature vector using a fine-tuned decoder.
[0019] Preferably, the first distance is Mahalanobis distance, and the second distance is Euclidean distance. Assigning non-uniform training weights based on the second distance between the current multidimensional drift feature vector and the multidimensional drift feature vector corresponding to the historical monitoring data includes: for each data segment contained in the historical monitoring data, extracting its corresponding multidimensional drift feature vector, calculating the Euclidean distance scalar value between the multidimensional drift feature vector and the current multidimensional drift feature vector; performing a negative exponential transformation on all calculated Euclidean distance scalar values and summing them to obtain a normalized proportional base; dividing the transformed value of each data segment by the normalized proportional base to generate a non-uniform training weight value assigned to each data segment in local update training.
[0020] Euclidean distance is used to measure the similarity between historical and current working conditions. A negative exponential transformation maps the distance to weight values, giving higher weights to historical samples with characteristics similar to the current working conditions, while the weights of samples with significant differences decay exponentially. After normalization, non-uniform training weights are generated for each historical data segment. This approach highlights historical experience closely related to the current drift state during local model updates, suppressing the interference of irrelevant old samples on parameter adjustments, thus achieving adaptive response to current geological changes while preserving historical knowledge.
[0021] Preferably, the step of adjusting the learning rate based on the first distance to perform a local update of the functional module to be updated includes: keeping the parameters in the initial control model, except for the functional module to be updated, absolutely frozen; feeding the augmented training set with the non-uniform training weights into the data generation model, evaluating the prediction residuals through the weighted mean square error loss function; executing the backpropagation algorithm to calculate only the network layer gradients of the functional module parameters to be updated, and using the weight decay optimizer to perform parameter update iterations on the network layer gradients based on the learning rate adjusted based on the first distance to complete the local update.
[0022] This invention provides an intelligent control system for mine grouting equipment, including a processor and a memory. The memory stores computer program instructions, and when the computer program instructions are executed by the processor, the aforementioned intelligent control method for mine grouting equipment is implemented.
[0023] By adopting the above technical solution, a computer program is generated from the above-mentioned intelligent control method for mine grouting equipment and stored in a memory so that it can be loaded and executed by a processor. In this way, a terminal device can be made based on the memory and the processor for convenient use.
[0024] The technical solution of the present invention has the following beneficial technical effects: By parallel computing the rate of change of higher-order statistics of monitoring data and the time-series gradient of the model feature SHAP attribution value, a multidimensional drift feature vector with high sensitivity to changes in working conditions is constructed. Mahalanobis distance is used to accurately determine distribution shifts. When the shift exceeds the limit, a drift pattern classifier identifies the specific geological cause and locates the functional module to be updated. Then, physical constraint equations conforming to the current fluid-solid correlation law are instantiated from a pre-set model library and integrated into virtual grouting samples driven by the conditional variational autoencoder generation mechanism. Subsequently, a negative exponential weighting strategy based on Euclidean distance is adopted to differentiate training weights for historical samples with different correlation degrees. The learning rate is dynamically adjusted according to the degree of drift, and only local parameter updates are performed on the affected modules. This method effectively solves the problems of traditional schemes such as insensitivity to early gradual drift detection, neglect of physical laws during the update process, model adaptation lag under small samples, and forgetting of historical effective knowledge caused by full retraining. It significantly improves the accuracy and timeliness of grouting control under complex mine geological conditions. Attached Figure Description
[0025] Figure 1 This is a flowchart of an intelligent control method for mine grouting equipment according to the present invention; Figure 2 This is a diagram illustrating the temporal gradient comparison of SHAP attribution values; Figure 3 This is a schematic diagram comparing the performance of ablation experiments. Detailed Implementation
[0026] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.
[0027] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.
[0028] This invention discloses an intelligent control method for mine grouting equipment, referring to... Figure 1 This includes the following steps: S1. Acquire monitoring data and generate multidimensional drift feature vectors.
[0029] Specifically, a multidimensional online monitoring data stream of the mine grouting process is acquired, and control commands for the grouting equipment are output by the initial control model. Based on the rate of change of higher-order statistics of the multidimensional online monitoring data stream within the sliding time window and the temporal gradient of the SHAP attribution value of the input features of the initial control model, a multidimensional drift feature vector is generated.
[0030] In the specific implementation process, a sensor network installed in the grouting pipeline and pump station collects time-series sensor data such as pressure, flow rate, and concentration as a multi-dimensional online monitoring data stream. The Apache Kafka stream processing platform is used to push the data to an in-memory database in real time. The PyTorch deep learning framework is used to build an initial control model, and the forward propagation algorithm is used to process the data stream, thereby calculating and outputting grouting equipment control commands to control the grouting pump speed and valve opening.
[0031] A fixed-size sliding time window is established in memory, progressing forward with each time step. The window length can be set according to the typical working condition change cycle of the mine grouting process, ranging from 300s to 600s; the time step size can be set according to the drift detection sensitivity requirements and computational resource conditions, ranging from 5s to 20s. The window slides forward continuously according to this time step, capturing time series segments of multi-dimensional online monitoring data streams such as grouting pressure, grouting flow rate, and grout density. In a single window containing 30 discrete samples, skewness and kurtosis calculation functions from a scientific computing library can be called to obtain the skewness and kurtosis values of each dimension of data within the current time window and the immediately preceding sliding time window. Then, a backward difference algorithm is used to subtract the corresponding previous skewness and kurtosis values from the current skewness and kurtosis values to obtain the difference. Taking grouting pressure as an example, if the current window skewness calculation value is 1.2 and the previous window storage value is 1, the difference between the two skewness values is 0.2. Dividing this difference by the time step of 10 seconds, we get the skewness change rate of 0.02 per second and the corresponding kurtosis change rate, which constitute the higher-order statistical change rate.
[0032] Meanwhile, for the initial control model constructed by the long short-term memory network, the deep interpreter algorithm can be used to perform inference-level black-box interpretation. A preset number of historical samples from stable operating cycles are selected as the background distribution basis. For each feature dimension of the input feature vector in the current time window, an independent baseline SHAP attribution value is calculated. For example, the baseline SHAP attribution value of the grouting hole pressure in the current time window is calculated to be 0.45, and the pipeline flow rate is 0.32.
[0033] Then, the historical baseline SHAP attribution value within the previous sliding time window is subtracted from the current attribution value. If the baseline value of pore pressure in the previous window is 0.4, then the first-order discrete difference of pore pressure is 0.05. The first-order discrete differences of all feature dimensions are divided by a time step of 10 seconds to obtain the corresponding temporal gradient values, such as the temporal gradient of pore pressure being 0.005 per second. These are then concatenated to generate the complete temporal gradient of the SHAP attribution value.
[0034] Finally, the rate of change of higher-order statistics and the temporal gradient of SHAP attribution values are linearly concatenated along the feature dimension to generate a multidimensional drift feature vector. For example... Figure 2 As shown in the figure, the importance of each input feature to the initial control model decision fluctuates. Among them, the time gradient of the attribution values of deep deformation and structural stress is significantly higher than that of other features, which are important feature sources driving model drift.
[0035] S2. Evaluate the features, distance, identification, drift causes, and modules.
[0036] Specifically, when the first distance between the multidimensional drift feature vector and the historical stable state feature distribution exceeds a preset threshold, the cause of the drift event and the functional modules to be updated in the initial control model are identified by the drift pattern classifier, and the appropriate constraint equations are obtained from the preset model library according to the cause.
[0037] In the specific implementation process, the first distance is represented by the Mahalanobis distance. The inverse matrix of the pre-calculated historical stable-state feature covariance matrix is read from the database. The Mahalanobis distance calculation function in the scientific computing library can be called to calculate the Mahalanobis distance between the current multidimensional drift feature vector and the historical stable-state feature mean vector. The preset threshold is based on the Mahalanobis distance distribution of the multidimensional drift feature vectors under historical stable states, and adopts... The principle is defined as follows: the threshold is the historical Mahalanobis distance mean plus three times the standard deviation. When this distance exceeds the preset threshold, the feature vector is input into the drift pattern classifier. In this embodiment, the Mahalanobis distance calculation function in the SciPy library is called, and the threshold is determined to be 2.58 through statistical analysis of 1000 hours of stable operating data.
[0038] The drift pattern classifier can be constructed using machine learning models such as multilayer perceptrons, support vector machines, or random forests, and is trained using a historical drift event dataset. Multidimensional drift feature vectors labeled with causes and modules to be updated are used as training samples, and supervised training is performed using the cross-entropy loss function. After training, the parameters are frozen for online inference. This classifier consists of an input layer, three fully connected hidden layers with 256, 128, and 64 nodes respectively, and a Softmax output layer with 6 nodes. After the feedforward input, it generates six predicted probability vector components for each category, such as outputs of 0.05, 0.12, 0.75, 0.03, 0.02, and 0.03. The Argmax function is used to extract the category index node containing the highest probability value of 0.75, which is then identified as label 2.
[0039] The pre-configured static dictionary hash table is constructed through the following steps: Collect samples of various historical drift events, label the corresponding flow field structure variation characteristics and the affected functional modules in the initial control model; establish a one-to-one mapping relationship between the category index of the drift mode classifier and the above-mentioned labeling information, and store it as a hash table structure. Use tag 2 as the key to perform a mapping query on the static dictionary hash table to obtain specific flow field structure variation characteristics as the identified causes, such as permeability step change caused by fractured rock mass penetration or local strong resistance in pipelines caused by sudden grout solidification. Simultaneously, obtain the parameter names of independent network layers by matching the field columns corresponding to the causes to locate the functional modules to be updated, such as the parameter matrix of the second long short-term memory hidden layer inside the initial control model.
[0040] Next, based on the identified causal strings, a traversal and matching process is performed in the offline configured pre-built model library to retrieve the general-level partial differential flow field parameter equations containing theoretical formulas such as Darcy's law from the exclusive model packages with matching names. Specifically: ; in, For permeability tensor; It is the fluid dynamic viscosity coefficient; Pressure distribution; This represents the grouting source flow rate. The grout viscosity sequence and instantaneous flow rate sequence are extracted from the current online monitoring data stream. Representative values are obtained after data smoothing, such as an average viscosity of 48 mPa. The instantaneous average flow rate is 150 L / min. Substituting these representative values into the corresponding fluid dynamic viscosity coefficient and dynamic flow boundary condition terms in the equations, we complete the constant-value assignment. After assignment, discretization and spatial topological mapping generate a set of constraint equations, which contain the following linear algebraic relationships: ; in, This is a vector of theoretical observations, which in this embodiment is a 2×1 column vector. Its function is to represent the theoretical observation results of the two sensors. This vector is obtained by multiplying the model state variables and the operator matrix. The observation operator matrix for dimension reduction mapping is used to achieve the mapping alignment from the spatial grid state to the actual position of the sensor. It is constructed by comparing and matching three-dimensional spatial coordinates. If the sensor coincides with the 500th grid, then the first row of this matrix is all 0 except for the 500th column which is 1. The state vector of all nodes in the grid is 100,000×1 and is used to represent the state distribution of the entire grid of the model. It is generated by the finite volume method.
[0041] The linear algebraic equation of the global network is: ;in, The global coefficient matrix is of size N×N and consists of the conductivity coefficients of each grid. It is used to constrain the neighborhood propagation relationship of the flow field equation. Its value is extracted and generated by the partial differential flow field equation mathematical discretization algorithm. This is a column vector, composed of grouting source terms, dynamic flow boundary conditions, and constant terms. It provides boundary constraints for solving the network node equations. Its values are obtained by real-time online monitoring of boundary data and filling with a replacement algorithm. The replacement algorithm specifically involves initializing the column vector. For a vector consisting entirely of zeros, the corresponding mesh node index is matched based on the three-dimensional spatial coordinates of the boundary grouting hole location, and the column vector is then... The element at that index position is replaced with a representative value of the instantaneous traffic volume from real-time online monitoring.
[0042] S3. Introduce mechanistic constraints to generate virtual grouting process data.
[0043] Specifically, after obtaining the appropriate constraint equations, the constraint equations are incorporated as loss terms into the data generation model to generate virtual grouting process data, and the corresponding virtual control commands are output by the initial control model.
[0044] In the specific implementation process, the data generation model employs a conditional variational autoencoder. The current online monitoring data stream is concatenated with a multidimensional drift feature vector and input into the encoder layer of this model, mapping and outputting the mean and variance vectors of a Gaussian latent variable distribution. Low-dimensional latent variable vectors are sampled from this data using reparameterization techniques, and their calculation relationship is as follows: ; in, The sampled low-dimensional latent variable vector serves as a hidden layer representation containing drift and variation features to pass data downstream. It is calculated from a combination of mean, standard deviation, and noise. The mean vector of the Gaussian latent variable distribution is used to represent the centroid of the latent space distribution and is obtained from the forward mapping output of the encoder network. The standard deviation vector of the Gaussian latent variable distribution is used to define the fluctuation range of the latent space distribution and is obtained by taking the square root of the variance vector output by the encoder. The sampling noise, which follows a standard normal distribution, serves to avoid the non-differentiability problem of the sampling operation during backpropagation. It is obtained by sampling through an environmental noise random generator. This is the symbol for element-wise multiplication between vectors.
[0045] The low-dimensional latent variable vector and the multi-dimensional drift feature vector are concatenated and input into the decoder layer. The result is then inversely projected to output the first batch of high-dimensional virtual full-space grid state tensors that contain the dimensions of all grid nodes. The virtual full-space grid state tensor Substituting directly into the aforementioned global network linear algebraic equilibrium equations, the root mean square error of the global physical equilibrium residuals is calculated. For example, if the calculated residual value is 0.12 MPa, it is dimensionless and multiplied by a preset constant, which is then set as the mechanism penalty function term. Then, combined with the aforementioned observation operator matrix The initial virtual monitoring data tensor aligned to the corresponding sensor positions is extracted from the high-dimensional virtual full-space grid state tensor using matrix multiplication, and used to calculate the subsequent reconstruction loss. Simultaneously, to constrain the latent variable distribution in the latent space of the data generation model, a KL divergence measure under the prior distribution is introduced, with the following relationship: ; in, For the KL divergence measure under the prior distribution, its function is to constrain the distribution of latent variables in the latent space to make it approach the standard normal distribution to avoid overfitting. It is approximated by internal summation and logarithmic integration. and These are the squares of the variance vector and the mean vector, respectively, both derived from the encoder's inference output.
[0046] ; in, The overall aggregated evaluation of the total loss serves as the sole objective function to guide the reverse update of the neural network weights, and is composed of three parts of loss. The native mean squared error reconstruction loss is used to measure the deviation between the decoder-generated tensor and the real online data. It is obtained by calculating the variance between the initial virtual monitoring data tensor and the real multidimensional online monitoring data stream. The preset constant coefficient, with a value ranging from 0.1 to 1, is used to adjust the weight ratio of the mechanism penalty function in the overall loss and is obtained by human experience in advance.
[0047] By reducing the total loss through backpropagation algorithm, the network is fine-tuned in multiple rounds. After reaching the preset 250-round cutoff period, the random environmental noise vector and multidimensional drift feature vector of standard normal distribution are extracted and input into the fine-tuned decoder. The decoder directly outputs virtual grouting process data that follows the characteristics of working condition variation and the physical mechanism of flow field. The data is converted and input into the initial control model to output virtual control commands such as adjusting the hydraulic pump gain to level 2.5.
[0048] S4. Integrate historical and virtual samples to construct an enhanced training set.
[0049] Specifically, an enhanced training set is constructed by integrating virtual grouting process data, corresponding virtual control commands, and pre-acquired historical monitoring data.
[0050] In the specific implementation process, the system reads historical monitoring data fragments stored in the time decay period. This data includes multi-dimensional feature sequences from the previous stable operation period and historical control commands issued. The time decay period can be set according to the rate of change of the mine's geological conditions, with a value range of 1 to 30 days. Historical data exceeding this period are not included in the training set due to excessive changes in geological conditions. The extracted historical monitoring data and the virtual grouting process data and their corresponding virtual control commands generated in step three are merged and concatenated using the concat function in the Pandas data processing library along the sample row dimension. This matrix concatenation operation can mix data from different periods and with different generation attributes to construct an enhanced training set that covers both baseline historical features and forward drift characteristics.
[0051] S5, Dynamically allocate weights and locally update the model output command.
[0052] Specifically, non-uniform training weights are assigned based on the second distance between the current multidimensional drift feature vector and the multidimensional drift feature vector corresponding to the historical monitoring data. The learning rate is adjusted in combination with the first distance to perform local updates of the functional modules to be updated. The updated control model outputs grouting equipment control commands.
[0053] In the specific implementation process, the second distance is manifested as Euclidean distance. For the offline historical samples in the augmented training set, the historical multidimensional drift feature vector corresponding to each segment is extracted in parallel. The spatial Euclidean norm difference between this vector and the current multidimensional drift feature vector is calculated to obtain the scalar value of the Euclidean distance. For example, the distance obtained for segment A, which is close to the current working condition, is 0.32, while the distance for segment B, which is far from the current condition, is 3.85. A transformation operation is performed using the Gaussian kernel decay formula, that is, the distance is multiplied by a hyperparameter factor set to 1.0 and then subjected to a negative exponential operation based on the natural logarithm base e. Through this transformation, the feature value of 0.32 is transformed into an index value of approximately 0.726, and 3.85 is transformed into approximately 0.021, thereby suppressing the influence of weakly correlated samples. The entire transformed sequence is arithmetically summed to obtain the normalized proportional base. Finally, the value of each segment is divided by this base, thereby generating and configuring non-uniform training weight values between 0 and 1. For virtually generated samples, the highest weight of one is assigned. Figure 3 As shown, the curve of this non-uniform training weight exhibits exponential decay, which can effectively cut off the historical knowledge forgetting caused by outdated heterogeneous data.
[0054] During the parameter iteration phase, an inversely proportional decay function between the Mahalanobis distance (reflecting drift intensity) and the baseline learning rate is established. The learning rate is adaptively adjusted in real-time before backpropagation using the LambdaLR scheduler in the PyTorch optimizer. Simultaneously, all parameters of other functional modules not identified in step two are absolutely frozen, and the weighted enhanced training set is fed into the initial control model in batches. After calculating the loss, the gradient of the network layers is obtained for the parameters of the functional modules to be updated identified in step two, and the AdamW weight decay optimizer is used for multiple iterations of updates according to the adjusted learning rate. After the local update is complete, the latest collected monitoring data flows through a fully connected mapping to output updated grouting pump operating frequency and pipeline valve opening control commands.
[0055] To verify the effects of introducing mechanistic constraints and local optimization, an experimental set containing 150,000 online monitoring data streams was constructed from a historical database. 4,500 time-series data streams induced by crack abrupt changes were extracted as the test set. Results showed that the virtual samples generated by the control group, relying solely on pure data-driven error, violated the law of mass conservation in 16.5% of cases, with a root mean square error (RMSE) of 0.52 MPa. This resulted in an updated action command accuracy of only 78.4% and a control delay lag of 14.2 s. In contrast, the experimental group employing the aforementioned physical mechanism constraints reduced the proportion of virtual violation samples to 1.2%, and the RMS error converged to 0.11 MPa. After local model updates, the accuracy of grouting action commands increased to 95.7%, and the control decision delay was reduced to 3.5 s. Figure 3The figure visually presents the performance comparison results of the two sets of experiments, verifying the effectiveness of the method of the present invention in dealing with sudden changes in complex working conditions.
[0056] This invention also discloses an intelligent control system for mine grouting equipment, including a processor and a memory. The memory stores computer program instructions, and when the computer program instructions are executed by the processor, an intelligent control method for mine grouting equipment according to the present invention is implemented.
[0057] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.
Claims
1. A method for intelligent control of mine grouting equipment, characterized in that, include: S1. Obtain multi-dimensional online monitoring data stream of the mine grouting process, and output grouting equipment control commands from the initial control model; Based on the rate of change of higher-order statistics of multidimensional online monitoring data stream within a sliding time window and the time-series gradient of SHAP attribution values of initial control model input features, a multidimensional drift feature vector is generated; S2, when the first distance between the multidimensional drift feature vector and the historical stable state feature distribution exceeds a preset threshold, the cause of the drift event and the functional modules to be updated in the initial control model are identified by the drift pattern classifier, and the appropriate constraint equations are obtained from the preset model library according to the cause; S3. Integrate the constraint equations as loss terms into the data generation model to generate virtual grouting process data. This includes: the data generation model being a conditional variational autoencoder; inputting the multidimensional online monitoring data stream and multidimensional drift feature vectors into the conditional variational autoencoder to generate the virtual full-space grid state tensor and the mean and variance of the latent variable distribution; subsequently, substituting the virtual full-space grid state tensor into the constraint equations to obtain the mechanism penalty function term, and extracting data from it through the observation operator matrix to calculate the original reconstruction loss; finally, adding the mechanism penalty function term, the original reconstruction loss, and the KL divergence term derived based on the mean and variance to obtain the total loss; updating the network parameters of the conditional variational autoencoder based on the total loss, and generating virtual grouting process data based on random noise and multidimensional drift feature vectors through a fine-tuned decoder; S4. The virtual control command is output from the initial control model; S5. The virtual grouting process data, the corresponding virtual control command, and the pre-acquired historical monitoring data are integrated to construct an enhanced training set; S6. Based on the second distance between the current multidimensional drift feature vector and the multidimensional drift feature vector corresponding to the historical monitoring data, non-uniform training weights are assigned, and the learning rate is adjusted in combination with the first distance to locally update the functional module to be updated. The updated control model outputs the grouting equipment control command.
2. The intelligent control method for mine grouting equipment according to claim 1, characterized in that, Based on the rate of change of higher-order statistics of the multidimensional online monitoring data stream within the sliding time window, the method includes: setting the time step and window length of the sliding time window, and acquiring the multidimensional online monitoring data stream of each stage in real time according to the time step; Calculate the skewness and kurtosis of each dimension in the multidimensional online monitoring data stream within the current sliding time window; use the backward difference algorithm to calculate the difference between the skewness and kurtosis of each corresponding dimension between the current sliding time window and the immediately preceding sliding time window; Divide each of the differences by the time step to obtain the skewness and kurtosis rates of the monitoring data for each dimension corresponding to each sliding time window, which are used as the rate of change of the higher-order statistics.
3. The intelligent control method for mine grouting equipment according to claim 1, characterized in that, The temporal gradient of the SHAP attribution value of the input features of the initial control model is obtained through the following steps: the inference stage of the initial control model is interpreted using an analytical algorithm, and the independent corresponding baseline SHAP attribution value is calculated for each feature dimension of the input features within the sliding time window. The first-order discrete difference is obtained by subtracting the benchmark SHAP attribution value of the corresponding feature dimension in the immediately preceding sliding time window from the benchmark SHAP attribution value of each feature dimension in the current sliding time window. Divide the first-order discrete differences of all feature dimensions by the time step of the sliding time window and form a vector, which is used as the temporal gradient of the SHAP attribution value of the input feature.
4. The intelligent control method for mine grouting equipment according to claim 1, characterized in that, The function module for identifying the cause of drift events and updating the initial control model through the drift pattern classifier includes: inputting the multidimensional drift feature vector into the drift pattern classifier and calculating the predicted probability vector components of each category label generated by the output layer. Extract the category index node containing the component with the largest numerical value; Read the dictionary hash table with pre-configured static association information, and query through the category index node mapping to obtain the specific flow field structure variation state characteristics as the cause of drift; By matching the field columns corresponding to the cause, the corresponding independent network layer parameter names are obtained to determine the functional module to be updated.
5. The intelligent control method for mine grouting equipment according to claim 1, characterized in that, The pre-set model library is a fluid-structure interaction model library. The step of identifying the causes of drift events and the functional modules to be updated in the initial control model through the drift pattern classifier, and obtaining the appropriate constraint equations from the pre-set model library according to the causes, includes: extracting partial differential flow field parameter equations with matching names and containing unassigned initial coefficient terms from the fluid-structure interaction model library according to the identified causes. Extract the slurry viscosity sequence and instantaneous flow rate sequence from the current multidimensional online monitoring data stream, and obtain the three-dimensional spatial coordinates corresponding to each grouting hole location and pipeline sensor; The temporal characteristics of the slurry viscosity sequence are substituted into the hydrodynamic viscosity coefficient term of the partial differential flow field parameter equation, and the instantaneous flow rate sequence is substituted into the dynamic flow rate boundary condition term of the corresponding spatial coordinate point. The parameter values are then assigned and the constraint equation is generated by discretization and spatial topology mapping.
6. The intelligent control method for mine grouting equipment according to claim 5, characterized in that, The process of generating the constraint equations through discretization and spatial topology mapping includes: calling a mathematical discretization algorithm component to transform the continuous partial differential flow field parameter equations with initial coefficients into a set of linear algebraic equilibrium equations between grid nodes according to the divided spatial grid. Based on the three-dimensional spatial coordinates of each grouting hole location and pipeline sensor, a dimension-reduced mapping observation operator matrix is constructed to form a set of constraint equations that match the current geological conditions.
7. The intelligent control method for mine grouting equipment according to claim 1, characterized in that, The first distance is Mahalanobis distance, and the second distance is Euclidean distance; non-uniform training weights are assigned based on the second distance between the current multidimensional drift feature vector and the multidimensional drift feature vector corresponding to the historical monitoring data, including: for each data segment contained in the historical monitoring data, extracting its corresponding multidimensional drift feature vector, and calculating the Euclidean distance scalar value between the multidimensional drift feature vector and the current multidimensional drift feature vector; Perform negative exponential transformation on all calculated Euclidean distance scalar values and sum them to obtain the normalized scaling factor; divide the transformed value of each data segment by the normalized scaling factor to generate a non-uniform training weight value assigned to each data segment in local update training.
8. The intelligent control method for mine grouting equipment according to claim 1, characterized in that, The step of adjusting the learning rate based on the first distance to perform local updates on the functional modules to be updated includes: keeping the parameters in the initial control model, except for the functional modules to be updated, absolutely frozen. The enhanced training set, which is assigned the non-uniform training weights, is fed into the data generation model, and the prediction residuals are evaluated using the weighted mean square error loss function. The backpropagation algorithm is executed to calculate only the network layer gradient of the functional module parameters to be updated. The weight decay optimizer is used to perform parameter update iteration on the network layer gradient based on the learning rate adjusted based on the first distance to complete the local update.
9. An intelligent control system for mine grouting equipment, characterized in that, include: A processor and a memory, wherein the memory stores computer program instructions that, when executed by the processor, implement the intelligent control method for mine grouting equipment according to any one of claims 1-8.
Citation Information
Patent Citations
Intelligent control device and control method for mine grouting equipment
CN112127918B
Coal mine dust diffusion simulation and control method based on artificial intelligence
CN119623238A
Gold mine vein measuring and positioning method based on multi-data coupling
CN121741896A