Microseismic full waveform reverse-time migration imaging roof safety warning and interlock control method

Abstract

The present application relates to the technical field of mining engineering and mine safety monitoring, and discloses a roof safety early warning and interlocking control method based on microseismic full waveform reverse-time migration imaging, which comprises the following steps: obtaining microseismic clusters through homology clustering, constructing a spatial anisotropic diffusion tensor and a sparse sensitivity matrix by using a previous period differential reflectivity tensor, combining with a tail wave inversion to obtain a velocity disturbance field and update a velocity model; performing reverse-time migration and wave field decoupling, calculating a current differential reflectivity tensor through cross-period wave field cancellation; weightedly constructing a roof damage evolution factor and its first-order and second-order change rates, outputting an early warning signal containing spatial positioning and triggering interlocking control when a composite accelerated evolution criterion is met; and feeding back the current tensor to the next period. The present application realizes dynamic and continuous tracking of roof damage and accelerated early warning of instability, and significantly improves the timeliness of mine safety response.

Description

Technical Field

[0001] This invention relates to the fields of mining engineering and mine safety monitoring technology, specifically to a method for roof safety early warning and interlocking control based on microseismic full-waveform reverse time migration imaging. Background Technology

[0002] During underground mining, the roof rock mass is susceptible to localized damage accumulation, fracture penetration, and sudden instability due to factors such as mining-induced stress disturbance, structural plane development, joint and fracture expansion, and stress redistribution. Roof instability not only affects the continuity of mining operations but can also cause equipment damage and personnel casualties. Therefore, how to predict the damage state of the stope roof in real time, accurately, spatially, and dynamically is a crucial technical issue in mine safety production.

[0003] Existing methods for monitoring the roof of mining areas mainly include stress monitoring, displacement monitoring, acoustic emission monitoring, microseismic location monitoring, and conventional geophysical imaging. Among these, microseismic monitoring, by collecting elastic wave signals released during rock fracturing, reflects the propagation of fractures and energy release within the rock mass and has been widely used for early warning of dynamic disasters in mines. However, traditional microseismic monitoring typically focuses on the location of microseismic events, the counting of events, energy release analysis, or source mechanism inversion, making it difficult to directly establish a quantitative relationship between microseismic waveform characteristics and the continuous damage state of the roof rock mass.

[0004] Furthermore, existing microseismic location or tomographic imaging methods are easily affected by background bedding, initial fractures, noise interference, and velocity model errors in complex mining environments, leading to a disconnect between imaging results and the actual damage evolution state. For roof damage, simply judging whether microseismic events have increased or energy has risen is often insufficient to identify whether the damage has entered an accelerated propagation stage, and it is also difficult to obtain the spatial distribution, evolution rate, and acceleration characteristics of local roof damage.

[0005] Reverse-time migration imaging (RTM) can image the subsurface medium structure using full waveform information and has a high wavefield reconstruction capability; wake interferometry is highly sensitive to minute velocity changes in the medium. However, current technologies typically use these two methods separately for imaging or velocity change analysis, lacking a cross-cycle closed-loop processing mechanism that feeds back imaging results from the previous monitoring cycle for inversion in the current cycle and further feeds back imaging results from the current cycle to the next cycle. Therefore, current technologies still suffer from problems such as difficulty in eliminating background bedding interference, disconnect between geophysical features and damage state, and difficulty in quantifying and identifying damage evolution rate and acceleration phases. The root cause lies in the lack of a cross-monitoring cycle closed-loop feedback mechanism, which isolates the inversion and imaging processes of each monitoring cycle, making it impossible to use historical damage identification results to guide the current cycle's processing to focus on areas that have changed. Summary of the Invention

[0006] The purpose of this invention is to provide a method for roof safety early warning and interlocking control based on microseismic full-waveform reverse time migration imaging, so as to solve the problems of existing microseismic monitoring methods that are difficult to dynamically quantify the spatial distribution and evolution speed of roof damage, and difficult to accurately identify the damage acceleration and instability stage and perform spatial positioning early warning.

[0007] This invention couples microseismic wake wave interferometry inversion, full-waveform reverse time migration imaging, differential reflectivity tensor, spatial anisotropic diffusion tensor, sparse sensitivity matrix, and damage evolution factor across cycles. It uses the imaging results of the previous cycle to constrain the velocity perturbation inversion and velocity model update of the current cycle, and feeds back the imaging results of the current cycle to the next cycle for processing, forming a closed-loop early warning and interlocking control mechanism. This improves the accuracy of identifying local damage to the roof, judging its evolution trend, and providing early warning of instability.

[0008] To achieve the above objectives, the present invention provides the following technical solution: a method for roof safety early warning and interlocking control based on microseismic full-waveform reverse time migration imaging, comprising the following steps:

[0009] S1: Discretize the target monitoring area into a spatial grid, configure the initial background wave velocity field, and set the monitoring cycle sequence;

[0010] S2: During the current monitoring period, the collected continuous waveform signals are dynamically truncated and clustered based on their similarity to obtain a cluster of microseismic events with similar origins.

[0011] S3: Obtain the differential reflectivity tensor calculated by reverse time-shift imaging in the previous monitoring cycle, construct a spatial anisotropic diffusion tensor by performing spectral decomposition on the differential reflectivity tensor, and reconstruct the spatial scattering sensitivity kernel function based on the diffusion tensor to establish a sparse sensitivity matrix characterizing the sensitivity relationship between wave velocity disturbance and wake phase delay.

[0012] S4: Using the wake sequence of the co-source microseismic cluster, combined with the sparse sensitivity matrix, inversion is performed to obtain the velocity relative disturbance field of the current monitoring period, and the three-dimensional velocity model is updated accordingly.

[0013] S5: Using the updated three-dimensional velocity model as input, perform reverse time migration on the microseismic waveforms of the co-source microseismic cluster, and decouple the elastic wave field during the reverse time migration process; based on the time difference cross-correlation imaging condition, perform cross-period wave field cancellation on the decoupled wave field of the current monitoring period and the previous monitoring period; wherein, the cross-period wave field cancellation includes matrix subtraction of the cross-correlation integral result of the current period decoupled wave field and the historical integral result of the previous monitoring period after event energy weighted standardization, so as to solve the differential reflectivity tensor of the current monitoring period;

[0014] S6: Based on the absolute value scalar of the velocity relative to the disturbance field and the Frobenius matrix norm of the differential reflectivity tensor, a multi-source joint weighted algebraic mapping is performed to construct the roof damage evolution factor; the first-order time rate gradient and the second-order time rate acceleration of the roof damage evolution factor are calculated; when the roof damage evolution factor and its first-order time rate gradient and second-order time rate acceleration satisfy the preset composite damage acceleration evolution criterion, an early warning signal of local roof instability acceleration containing spatial positioning code is output, and a hardware-level interrupt or interlocking control action is triggered through the mine safety control bus. The interlocking control action includes stopping the operation of mining equipment or cutting off the power supply to the designated working area.

[0015] Furthermore, the differential reflectivity tensor of the current monitoring period is stored as the input tensor for constructing the spatial anisotropic diffusion tensor in step S3 of the next monitoring period, so as to form a feedback loop across monitoring periods.

[0016] Furthermore, when the current monitoring period is the initial monitoring period, the differential reflectance tensor of the previous monitoring period is initialized to a preset three-dimensional all-zero tensor, or initialized to a reference reflectance tensor obtained based on the complete roof rock mass background model.

[0017] Further, in step S3, a spatially anisotropic diffusion tensor is constructed by performing spectral decomposition on the differential reflectivity tensor, specifically as follows:

[0018] The differential reflectivity tensor, which is a three-dimensional second-order symmetric tensor, is subjected to eigenvalue decomposition to obtain eigenvalue scalars and corresponding unit eigenvectors;

[0019] Construct the corresponding directional projection tensor based on the outer product of each unit eigenvector and its transpose, and perform a weighted summation of each directional projection tensor according to the corresponding eigenvalues. Then, perform an algebraic summation of the summation result with the three-dimensional unit matrix scaled by the preset background isotropic diffusion constant to obtain the spatial anisotropic diffusion tensor.

[0020] Furthermore, the event energy weighted normalization process includes:

[0021] If there are no historical events corresponding to the current microseismic cluster in the previous monitoring period, then the all-zero tensor is used as the historical integration result.

[0022] If multiple historical events corresponding to the current co-source microseismic cluster exist within the previous monitoring period, then the historical conjugate wavefield tensor integral results of each event are weighted and averaged using the scalar seismic moments of each event as weights to obtain the standardized historical integral results.

[0023] Furthermore, the top plate damage evolution factor mentioned in step S6 satisfies:

[0024] ;

[0025] in, The roof damage evolution factor; The relative perturbation value of wave velocity in the velocity relative perturbation field in the grid element The absolute value scalar at that point; The differential reflectivity tensor in the mesh element The Frobenius matrix norm at that location; and These are the preset weighting coefficients.

[0026] Furthermore, the preset weighting coefficients and The data were obtained through uniaxial graded loading damage evolution tests on homogeneous rock specimens combined with least squares fitting calibration, and then corrected by scale mapping through on-site calibration tests in the mining area.

[0027] Furthermore, the first-order time-rate gradient and the second-order time-rate acceleration mentioned in step S6 are calculated using the following discrete difference formulas:

[0028] ,

[0029] ;

[0030] in, The preset monitoring cycle step size.

[0031] Furthermore, the criterion for accelerated evolution of composite damage is:

[0032] ;

[0033] in, To preset the ultimate rock mass fracture damage threshold, the The critical damage value for macroscopic fracture penetration was determined based on uniaxial graded loading damage evolution tests on homogeneous rock specimens.

[0034] Furthermore, the updating of the three-dimensional velocity model in step S4 specifically involves: recursively updating the three-dimensional velocity model of the previous monitoring period based on the velocity relative disturbance field, so that the three-dimensional velocity model continuously evolves in the monitoring period sequence.

[0035] Furthermore, the inversion in step S4 is solved using regularization constraints, which include at least one of spatial smoothing constraints, Tikhonov regularization constraints, or anisotropy constraints based on geological structure orientation.

[0036] Compared with the prior art, the present invention has the following beneficial effects:

[0037] First, this invention utilizes the differential reflectivity tensor obtained through reverse time migration imaging in the previous monitoring cycle to construct a spatial anisotropic diffusion tensor, and reconstructs the spatial scattering sensitivity kernel function based on the spatial anisotropic diffusion tensor. This allows the current cycle's wake interferometry inversion to be spatially constrained by the physical information of the crack structure changes in the previous cycle, automatically focusing on the identified damage area and its dominant direction, thereby improving the sensitivity of the velocity relative disturbance field inversion results to the roof damage structure.

[0038] Second, the present invention updates the three-dimensional velocity model with the velocity relative to the disturbance field obtained by the inversion of the current period, and further uses it for the full waveform reverse time migration imaging of microseismic events in the current period, so that the velocity model update is coupled with the imaging process, which can reduce the influence of fixed background bedding and velocity model errors on the imaging results.

[0039] Third, by using time-difference cross-correlation imaging conditions and cross-period wavefield cancellation, this invention performs matrix subtraction between the cross-correlation integral result of the decoupled wavefield in the current period and the standardized historical integral result of the previous monitoring period to obtain the differential reflectivity tensor of the current monitoring period, thereby highlighting the newly added or enhanced crack reflection, scattering and damage areas in the current period.

[0040] Fourth, this invention performs multi-source joint weighted algebraic mapping between the absolute value scalar of velocity relative to the perturbation field and the Frobenius matrix norm of the differential reflectivity tensor to construct a roof damage evolution factor. This allows wave velocity variation information and reflection structure variation information to participate in damage characterization from two physical dimensions: "medium property variation" and "structural interface variation," thereby improving the reliability and noise resistance of roof damage quantification results.

[0041] Fifth, this invention not only determines whether the roof damage evolution factor exceeds the threshold, but also further calculates the first-order time rate gradient and the second-order time rate acceleration, and uses a composite damage acceleration evolution criterion to identify whether local roof damage has entered the accelerated expansion stage. It can effectively distinguish between stable high damage state and accelerated instability state, thereby significantly reducing the false alarm rate caused by the natural fluctuations of local roof damage, and gaining a valuable early warning window for mine safety evacuation and emergency response.

[0042] Sixth, this invention stores the differential reflectivity tensor of the current monitoring cycle as the input tensor for constructing the spatial anisotropic diffusion tensor of the next monitoring cycle, forming a differential reflectivity tensor feedback closed-loop mechanism across monitoring cycles, thereby realizing continuous tracking and dynamic prediction of the entire process of damage evolution of the mining roof from initiation to expansion and accelerated instability. Attached Figure Description

[0043] Figure 1This is a diagram illustrating the method steps of the present invention.

[0044] Figure 2 This is a schematic diagram of the data flow in the cross-cycle feedback closed loop of this invention.

[0045] Figure 3 This is a logic diagram of hierarchical early warning interlocking control in an embodiment of the present invention. Detailed Implementation

[0046] The technical solution of the present invention will be further described in detail below with reference to specific embodiments. It should be understood that the following embodiments are only used to explain the present invention and do not constitute a limitation on the scope of protection of the present invention. Based on the technical concept of the present invention, any equivalent substitutions or adaptive adjustments made by those skilled in the art to the order of steps, parameter selection, algorithm model, sensor type, server architecture, or early warning control method should fall within the scope of protection of the present invention.

[0047] Example 1

[0048] This embodiment provides a method for roof safety early warning and interlocking control based on microseismic full-waveform reverse time migration imaging. This method can be applied to the roof of underground mines, surrounding rock of roadways, boundary rock mass of goaf areas, or other rock mass engineering areas that require dynamic damage monitoring and safety early warning.

[0049] In this embodiment, continuous microseismic waveform signals are collected by a distributed network of microseismic detectors deployed around the target monitoring area. The signals are then processed by edge computing nodes or a central server, including waveform truncation, event identification, homology clustering, kernel function dynamic reconstruction, wake wave interferometry inversion, reverse time migration imaging, cross-period wave field cancellation, differential reflectivity tensor solution, damage evolution factor construction, and early warning criterion output.

[0050] The method specifically includes the following steps.

[0051] S1: Discretize the target monitoring area into a spatial grid, configure the initial background wave velocity field, and set the monitoring cycle sequence.

[0052] First, a three-dimensional spatial mesh model is established for the target monitoring area. This three-dimensional spatial mesh model consists of multiple mesh elements, each corresponding to a spatial location within the target monitoring area.

[0053] In one implementation, the three-dimensional spatial mesh model can be divided using a regular hexahedral mesh; in other implementations, irregular tetrahedral meshes or hybrid meshes can be used, depending on the morphology of the stope boundary, the contour of the roadway, the distribution of fault structures, or the occurrence of ore bodies. The selection of the mesh scale needs to comprehensively consider the spacing of the microseismic detectors, the target monitoring accuracy, and the computing resources of the central server. For example, the mesh scale can be set to 1m to 5m, with a finer density of 0.5m for key monitoring areas and a wider density of 10m for peripheral areas.

[0054] Preferably, the spatial grid step size should satisfy the grid dispersion constraint condition for numerical wave field propagation, that is, the spatial step size should not be greater than the ratio of the minimum propagation velocity of elastic waves in the target area to twice the highest effective cutoff frequency, so as to reduce the grid dispersion error in the extrapolation process of finite difference wave field.

[0055] After discretizing the spatial grid, the initial background wave velocity field of the target monitoring area is configured. The initial background wave velocity field can be obtained through methods including, but not limited to: field acoustic logging data, borehole wave velocity test data, historical microseismic location inversion results, geological exploration profile data, or laboratory ultrasonic rock mechanics test results. The initial background wave velocity field may include P-wave velocity parameters and S-wave velocity parameters at each grid cell, or it may only include P-wave velocity parameters, with the S-wave velocity parameters calculated based on a preset empirical relationship between P-wave and S-wave velocity ratios.

[0056] At the same time, a monitoring cycle sequence is set:

[0057] ;

[0058] in, Indicates the first At the time corresponding to each monitoring cycle, the adjacent monitoring cycles satisfy the following conditions: , The monitoring cycle step is preset. The setting of the monitoring cycle step needs to comprehensively consider the mining advance speed, the frequency of microseismic activity, the roof stability level, and the real-time computing capabilities of the central server. For example, It can be set to 5 minutes, 10 minutes, 15 minutes, 30 minutes or 1 hour. In areas with high stress concentration or strong mining disturbance, a shorter monitoring cycle should be used.

[0059] S2: During the current monitoring period, the collected continuous waveform signals are dynamically truncated and clustered based on their similarity to obtain a cluster of microseismic events with similar origins.

[0060] During the current monitoring period, the distributed microseismic geophone network continuously acquires continuous waveform signals within the target monitoring area. The distributed microseismic geophone network may include multiple three-component microseismic geophones, single-component geophones, or distributed fiber optic acoustic sensing units. Each geophone establishes a communication connection with the central server via wired Ethernet, Wi-Fi wireless network, or downhole industrial ring network.

[0061] After receiving continuous waveform signals, the central server or edge computing nodes first preprocess them. Preprocessing operations include, but are not limited to: removing the DC component and linear trend term of the signal; bandpass filtering (the passband range is set, for example, from 10Hz to 500Hz, specifically determined based on the dominant frequency range of the microseismic signal in the target monitoring area); automatic identification and removal of abnormal channels; amplitude normalization of each channel; and multi-channel time synchronization correction based on GPS or NTP protocols. The purpose of preprocessing is to improve the sensitivity of subsequent event-triggered detection and the reliability of waveform similarity judgment.

[0062] Subsequently, dynamic event capture is performed on the preprocessed continuous waveform signal. In one embodiment, a long-short time average ratio operator is used to perform a streaming sliding window scan on the waveform data of each channel. The short time window length can be set to 0.1s to 1s, the long time window length can be set to 1s to 10s, and the long-short time window ratio threshold can be set to 2 to 5. When the number of triggering channels is greater than the preset network triggering threshold (e.g., at least 4 channels triggering simultaneously), a microseismic event is confirmed. The initial location of the triggering event can be determined by using a travel-time grid search method based on finite difference decomposition of the process function equation, or an absolute location method based on ray tracing, to determine the grid cell to which the initial seismic source location belongs.

[0063] For the detected microseismic events, the corresponding displacement timing signals are captured at each detector channel. The start time of the capture window is set to a preset buffer period (e.g., 0.05s to 0.2s) before the theoretical arrival time of the P-wave. The end time of the capture window is determined based on the attenuation of the root mean square energy envelope of the signal moving window to twice the background noise level, ensuring that the captured segment simultaneously includes the first arrival wave, the main phase, and a sufficiently long wake wave signal.

[0064] After obtaining multiple microseismic waveform segments, the segments are clustered based on their similarity to the source. In this invention, the similarity clustering refers to the process of grouping microseismic events into clusters based on waveform similarity and spatial proximity. Specifically, it includes: for any two microseismic events detected within the current monitoring period... and Calculate the normalized cross-correlation function of each cross-correlation function on each normally operating detector channel, take the maximum absolute value of the cross-correlation coefficient of each channel on the time delay axis, and calculate its arithmetic mean or median to obtain the average maximum cross-correlation coefficient. Simultaneously calculate the initial spatial distance between the two events. When both conditions are met and When, determine the event and Satisfying the homology condition, among which To preset the similarity threshold, This is a preset distance threshold.

[0065] Based on the aforementioned homology assessment results, density clustering, hierarchical clustering, or spectral clustering algorithms can be used to cluster all microseismic events detected within the current monitoring period, outputting several homogeneous microseismic clusters. Each homogeneous microseismic cluster contains at least two microseismic events. These clusters serve as the basic data processing units for subsequent wake-wave interferometric inversion and reverse-time migration imaging. Because events within homogeneous microseismic clusters exhibit high waveform similarity and comparable propagation paths, the reliability of wake-wave phase delay extraction and the stability of velocity relative perturbation inversion can be significantly improved.

[0066] S3: Obtain the differential reflectivity tensor calculated by reverse time migration imaging in the previous monitoring cycle, construct a spatial anisotropic diffusion tensor by performing spectral decomposition on the differential reflectivity tensor, and reconstruct the spatial scattering sensitivity kernel function based on the diffusion tensor to establish a sparse sensitivity matrix characterizing the sensitivity relationship between wave velocity disturbance and wake phase delay.

[0067] In the current monitoring cycle Initially, the central server reads the data from the storage unit for the previous monitoring period. Differential reflectivity tensor calculated and stored using reverse time-lapse imaging This differential reflectivity tensor is used to characterize the spatial distribution characteristics of reflective interface changes, crack propagation, abrupt changes in wave impedance, or local scattering anomalies in the target monitoring area during the previous monitoring period.

[0068] In this invention, the differential reflectivity tensor of the previous monitoring cycle is not an isolated imaging result, but rather serves as input data for the dynamic reconstruction of the kernel function in the current monitoring cycle. That is, the imaging differential result of the previous cycle is used to physically adjust the spatial distribution of the spatial scattering sensitivity kernel function in the current cycle, so that the wake inversion in the current cycle can automatically focus on the fracture development region and its dominant direction already identified in the previous cycle.

[0069] When the system is in the initial monitoring cycle of its first operation, since the differential reflectivity tensor calculated by reverse time migration imaging in the previous monitoring cycle does not yet exist, the differential reflectivity tensor of the previous monitoring cycle can be initialized to a preset three-dimensional all-zero tensor; or, based on the complete roof rock mass background model in an undisturbed state, a reference reflectivity tensor is obtained through a single reference reverse time migration imaging, and this reference reflectivity tensor is used as the input tensor for constructing the spatial anisotropic diffusion tensor in the initial monitoring cycle.

[0070] To obtain the differential reflectance tensor of the previous monitoring period Subsequently, spectral decomposition is performed on it to construct a spatially anisotropic diffusion tensor. Specifically, the differential reflectivity tensor is a three-dimensional second-order symmetric tensor, and its eigenvalue decomposition is performed on it:

[0071] ;

[0072] in, Indicates the first A scalar of eigenvalues; This represents the unit eigenvector corresponding to the eigenvalue; The directional projection tensor, constructed from the outer product of a unit eigenvector and its transpose, describes the damage structure in the th... Spatial orientation information in each main direction.

[0073] Based on the above eigenvalue decomposition results, the corresponding directional projection tensors are constructed using the outer product of each unit eigenvector and its transpose. Then, each directional projection tensor is weighted according to its own characteristics. Perform a weighted summation, where To preset the damage sensitivity adjustment coefficient, the summation result is then compared with the preset background isotropic diffusion constant. Scaled 3D identity matrix By performing algebraic summation, we obtain the spatial anisotropic diffusion tensor for the current monitoring period:

[0074] ;

[0075] in, This represents the spatial anisotropic diffusion tensor constructed during the current monitoring period; The preset background isotropic diffusion constant is used to characterize the baseline diffusion capacity of undamaged rock mass; This is the damage sensitivity adjustment coefficient. This function form makes the diffusion enhancement of the spatial anisotropic diffusion tensor in the main damage direction proportional to the magnitude of the differential reflectivity change in that direction, thereby effectively converting the intensity information of the damage area identified in the previous cycle into the focusing weight of the sensitivity kernel in the current cycle. The value can be calibrated according to the rock mass type and monitoring accuracy requirements. An example value is 500 to 2000.

[0076] The above construction method enables the spatial anisotropic diffusion tensor The enhanced diffusion coefficient in the main direction of the damage structure in the previous cycle guides the wake scattering sensitivity kernel function of the current cycle to physically focus on the known damage region and its dominant direction.

[0077] After obtaining the spatially anisotropic diffusion tensor, the spatial scattering sensitivity kernel function is reconstructed based on this diffusion tensor. Spatial scattering sensitivity kernel function is used to quantitatively describe different grid elements. Wave velocity disturbance at the observation point The degree of contribution of the tail wave phase delay. In one implementation, the spatial anisotropic diffusion tensor is used. Substituting the medium property matrix into the control equations for the anisotropic transport of elastic wave energy, numerical solutions and spatiotemporal convolution integration are performed. Specifically, spatiotemporal convolution integration is performed on the forward energy density field propagating from the microseismic source location to any grid element, and the accompanying energy density field propagating backward from the receiver location to any grid element. The integration results are then normalized using the total wave field energy density scalar propagating from the receiver location to the microseismic source location, thereby establishing the spatial coordinates of the grid element. Spatial sensitivity mapping of phase delay of multiple scattering tail waves.

[0078] Furthermore, based on the discrete values ​​of the spatial scattering sensitivity kernel function at each grid cell, the discrete coefficients of the sensitivity kernel function are matrix-assembled to establish a sparse sensitivity matrix characterizing the sensitivity relationship between wave velocity perturbation and wake phase delay. The rows of this sparse sensitivity matrix correspond to different pairs of microseismic events, station combinations, or wake time windows, while the columns correspond to different grid elements within the target monitoring area. The matrix elements... The numerical representation of the first The relative perturbation of wave velocity in the first grid cell affects the first... The contribution of each wake phase delay observation. For matrix elements smaller than a preset sparsity threshold. The portion can be directly set to zero to save storage and computational costs.

[0079] S4: Using the wake sequence of the co-source microseismic cluster, combined with the sparse sensitivity matrix, inversion is performed to obtain the velocity relative disturbance field of the current monitoring period, and the three-dimensional velocity model is updated accordingly.

[0080] Using the wake sequences of each microseismic event in the same microseismic cluster obtained in step S2, combined with the sparse sensitivity matrix established in step S3 Regularized inversion is performed to solve for the velocity relative perturbation field of the current monitoring period.

[0081] First, the wake phase delay observation is extracted. For microseismic event pairs within a co-originating microseismic cluster, waveform segments within their wake time windows are selected. The start time of the wake time window should ensure that it avoids the direct wave and the main reflected phase, typically taken after 2 to 5 times the S-wave travel time from the time of origin; the end time of the wake time window can be dynamically determined based on the signal-to-noise ratio. In one implementation, the relative velocity change of the wake sequence of the co-originating microseismic event pair is extracted using the time-axis stretching method: by applying a time-axis stretching fine-tuning operator to the wake signal of one of the events for resampling, its cross-correlation coefficient with the wake signal of the reference event is maximized, and the relative velocity change observation scalar is obtained by converting it with the optimal stretching factor. The relative velocity change observation scalars corresponding to multiple event pairs, multiple station channels, and multiple wake time windows are concatenated and assembled into a known observation vector. .

[0082] Establish the linearized inversion equation:

[0083] ;

[0084] in, This is the observation vector for the tailwave phase delay; The sparse sensitivity matrix established for step S3; Let be the velocity-relative perturbation field vector to be solved; This represents the observation error term.

[0085] Since inversion problems are often underdetermined or mixed-determined, regularization constraints are needed to improve the stability and uniqueness of the inversion. The following regularized least squares objective function is constructed:

[0086] ;

[0087] The first term is the data fitting term; the second term is the regularization constraint term. For spatial smoothing operators (such as first- or second-order spatial difference operators). The first term is the space smoothing regularization parameter; the second term is the cross-period stability constraint term. These are cross-cycle stability constraint parameters used to suppress non-physical abrupt changes in the velocity disturbance field between adjacent monitoring cycles. Regularization parameter. and It can be adaptively determined using the L-curve method or the generalized cross-validation method.

[0088] In other implementations, the regularization term Tikhonov regularization constraints, total variational sparsity constraints, or anisotropic smoothing constraints based on geological structure orientation can also be used in replacement or combination.

[0089] The objective function described above can be efficiently solved using the conjugate gradient algorithm, the LSQR algorithm, or other iterative methods suitable for large sparse linear equation systems, reconstructing the relative velocity perturbation field at each grid cell within the target monitoring area. .

[0090] velocity relative to the perturbation field It characterizes the minute relative change in elastic wave velocity at each grid cell within the current monitoring period compared to the wave velocity in the previous monitoring period. Its spatial distribution indirectly reflects the changes in local medium properties within the roof rock mass caused by damage evolution, fracture propagation, or stress adjustment.

[0091] The three-dimensional velocity model is updated based on the velocity relative to the perturbation field. In a preferred embodiment, a recursive update strategy is employed, using the three-dimensional velocity model from the previous monitoring period... As an update benchmark:

[0092] ;

[0093] in, This is the updated 3D velocity model for the current monitoring cycle; This is the velocity perturbation mapping coefficient, typically set to 1. The advantage of using a recursive update method is that the velocity model can continuously evolve over the monitoring cycle sequence, tracking in real time the wave velocity change trend caused by progressive damage to the roof rock mass.

[0094] For cases where both P-wave and S-wave propagation are considered simultaneously, P-wave velocity models and S-wave velocity models can be established separately and then inverted and updated separately; or only the P-wave velocity relative disturbance field can be inverted, and the update amount of the S-wave velocity model can be calculated based on the preset dynamic P-wave to S-wave velocity ratio relationship.

[0095] The updated three-dimensional velocity model will be used as the input wave velocity model for reverse time migration imaging in step S5, so that the reverse time migration can adapt to the actual wave velocity distribution corresponding to the current damage state of the roof rock mass, thereby effectively improving the imaging positioning accuracy and the ability to identify reflection anomalies.

[0096] S5: Using the updated three-dimensional velocity model as input, perform reverse time migration on the microseismic waveforms of the co-source microseismic cluster, and decouple the elastic wave field during the reverse time migration process; based on the time difference cross-correlation imaging condition, perform cross-period wave field cancellation on the decoupled wave field of the current monitoring period and the previous monitoring period; wherein, the cross-period wave field cancellation includes matrix subtraction of the cross-correlation integral result of the current period decoupled wave field and the standardized historical integral result of the previous monitoring period to calculate the differential reflectivity tensor of the current monitoring period.

[0097] The updated three-dimensional velocity model obtained in step S4 As the input wave velocity model, full-waveform reverse time migration imaging is performed on the microseismic waveform data of the same source microseismic cluster.

[0098] During the reverse-time migration, wavefield extrapolation is performed based on the three-dimensional two-way elastic first-order velocity-stress wave equations. Specifically, the displacement waveform time series of newly occurring microseismic events belonging to the same source microseismic cluster acquired within the current monitoring period are used as source boundary conditions or detector boundary conditions, and the source wavefield is extended forward along the time axis. Simultaneously, the waveforms recorded by each detector are time-reversed and used as boundary conditions, and the detector wavefield is extended backward along the time axis. Wavefield extrapolation can be numerically implemented using the higher-order finite difference method, pseudospectral method, or finite element method.

[0099] Within each time step of the wavefield extrapolation propagation loop, an elastic wavefield decoupling operation is performed. The purpose of elastic wavefield decoupling is to separate the coupled elastic vector wavefield into longitudinal (P-wave) and transverse (S-wave) components, thereby reducing imaging crosstalk between different wave types and improving the imaging resolution of damaged structures on the top plate (such as cracks and fissures). In one implementation, the spatial divergence differential operator is utilized. Acting on an elastic vector wave field, the scalar potential field of the longitudinal wave is extracted; the spatial curl differential operator is used. The transverse wave vector potential is extracted by applying it to the elastic vector wave field. Then, the spatial gradient operator of the potential field is used... By using the spatial inverse curl operator, the elastic longitudinal wave vector wave field components and the elastic transverse wave vector wave field components are reconstructed.

[0100] Time-difference cross-correlation imaging conditions are applied to the decoupled wavefield components. Let the current monitoring period be... The decoupled forward extension wave field is The decoupled reverse extension wave field is The cross-correlation integral result of the current periodic decoupled wavefield is:

[0101] ;

[0102] in, Record the duration of the waveform; The result of the cross-correlation integral for the current monitoring period is essentially the original reflectivity imaging tensor without cross-period cancellation.

[0103] Before performing cross-period wavefield cancellation, the historical integration results from the previous monitoring period need to be standardized. The purpose of standardization is to make historical imaging results from different monitoring periods comparable on an energy scale, thereby enabling cross-period difference to effectively remove the wavefield response of constant geological background bedding.

[0104] The specific rules for standardization are as follows:

[0105] (1) If there are no historical events corresponding to the current microseismic cluster in the previous monitoring period (i.e., the cluster is appearing for the first time), then the all-zero tensor is used as the historical integration result. At this point, the cross-period wavefield phase degenerates into the direct output of the current period imaging result, i.e. .

[0106] (2) If multiple historical events corresponding to the current co-source microseismic cluster exist within the previous monitoring period, the scalar seismic moments of each historical event are used as weights to perform a weighted average of the historical conjugate wavefield tensor integral results of each event, resulting in the standardized historical integral results:

[0107] ;

[0108] in, The results are the standardized historical integrals. The first in the previous monitoring period The result of the conjugate wavefield tensor integral of a historical event; For the first Scalar seismic moments of historical events; This represents the total number of historical events corresponding to the same microseismic cluster within the previous monitoring period.

[0109] The advantage of using scalar seismic moments as weights for weighted averaging is that it can highlight the contribution weight of high-energy events to the background medium structure, making the standardized historical integral results more representative of the main geological background bedding response characteristics corresponding to the cluster, thereby more effectively stripping away stable background structure information during subsequent cross-period cancellation.

[0110] Subsequently, the cross-correlation integral result of the decoupled wavefield in the current period is subtracted from the standardized historical integral result of the previous monitoring period by matrix subtraction to calculate the differential reflectivity tensor of the current monitoring period:

[0111] ;

[0112] in, This is the differential reflectance tensor for the current monitoring period.

[0113] The aforementioned cross-period wavefield destructive processing effectively suppresses the contribution of constant background structures such as stable geological stratification interfaces and unchanged regions in the mining roof to the imaging results, while significantly enhancing newly added or enhanced fracture propagation regions, local scattering anomalies, and reflectivity change regions within the current monitoring period compared to the previous monitoring period. This provides purer and more sensitive structural response information for constructing the roof damage evolution factor in subsequent step S6.

[0114] S6: Based on the absolute value scalar of the velocity relative to the disturbance field and the Frobenius matrix norm of the differential reflectivity tensor, a multi-source joint weighted algebraic mapping is performed to construct the roof damage evolution factor; the first-order time rate gradient and second-order time rate acceleration of the roof damage evolution factor are calculated; when the roof damage evolution factor and its first-order time rate gradient and second-order time rate acceleration satisfy the preset composite damage acceleration evolution criterion, an early warning signal of accelerated local instability of the roof containing spatial positioning code is output, and a hardware-level interrupt or interlocking control action is triggered through the mine safety control bus; and the differential reflectivity tensor of the current monitoring cycle is stored as the input tensor used to construct the spatial anisotropic diffusion tensor in step S3 of the next monitoring cycle, so as to form a feedback closed loop across monitoring cycles.

[0115] The velocity relative perturbation field for the current monitoring period is obtained. and differential reflectivity tensor Subsequently, multi-source joint weighted algebraic mapping is performed to construct the roof damage evolution factor.

[0116] For each grid cell element Top plate damage evolution factor Defined as:

[0117] ;

[0118] in, For the velocity relative perturbation field in the grid element The absolute value scalar at that location reflects the relative variation of the rock mass wave velocity at that location, and indirectly characterizes the degree of medium damage caused by the initiation, expansion, or change in porosity of the rock mass. For differential reflectivity tensor The Frobenius matrix norm is defined as follows: Its physical meaning is the intensity of reflected or scattered energy that is increased or enhanced in the current cycle compared to the previous cycle; and These are preset weighting coefficients used to adjust the relative contribution weights of wave velocity disturbance information and reflection structure change information in damage characterization.

[0119] The technical advantages of the aforementioned multi-source joint weighted algebraic mapping are as follows: the velocity relative to the perturbation field is more sensitive to the overall change in the density of microfractures within the rock mass, enabling it to capture the gradual evolution of damage at the medium scale; while the Frobenius norm of the differential reflectivity tensor is more sensitive to the abrupt changes in reflectivity caused by the formation and expansion of local fracture surfaces, enabling it to capture sudden changes in damage at the structural scale. These two aspects complement and synergistically constrain each other, allowing the roof damage evolution factor to comprehensively characterize the damage state of the roof rock mass from two physical dimensions: "changes in medium properties" and "changes in structural interfaces," significantly improving the reliability and noise resistance of the damage quantification results.

[0120] Furthermore, in order to identify whether the damage has entered the accelerated propagation stage, the first-order time rate gradient of the roof damage evolution factor on the monitoring period sequence was calculated. and second-order time rate of change acceleration Both are calculated using the following discrete difference formulas:

[0121] ;

[0122] ;

[0123] in, The preset monitoring cycle step size. Characterizes the rate of damage evolution, The acceleration characterizes the rate of damage evolution, i.e. whether the damage is in a state of accelerated growth.

[0124] Based on the above indicators, the following composite damage acceleration evolution criterion is used to determine the early warning of local instability of the roof. When the same grid element In the current monitoring cycle When all three of the following conditions are met, the region corresponding to the grid cell is determined to have entered a state of accelerated damage evolution:

[0125] ;

[0126] in, The preset ultimate rock mass fracture damage threshold is defined as follows: The determination method is as follows: a uniaxial graded loading damage evolution test is carried out on a rock specimen with the same lithology as the roof of the target monitoring area. The critical value of the roof damage evolution factor when the specimen enters the macroscopic crack penetration stage is recorded, and then multiplied by the preset safety reduction coefficient.

[0127] The physical meaning and technical effect of this composite criterion are as follows: First condition Ensure that the cumulative damage has reached a critical level to eliminate false alarms caused by minor damage fluctuations; second condition Ensure the damage is in a continuous growth phase rather than a recovery or stabilization phase; third condition This invention specifically identifies the physical critical point at which damage evolution transitions from uniform development to accelerated instability. Compared to existing technologies that rely solely on a single threshold of absolute damage value, the composite criterion of this invention effectively distinguishes between a "stable high-damage state" (high damage value but no longer increasing or increasing at a deceleration) and an "accelerated instability state" (high damage value and accelerating growth). This significantly reduces the false alarm rate caused by natural fluctuations in local roof damage, providing a valuable early warning window for mine safety evacuation and emergency response.

[0128] When a grid cell is determined to have entered a state of accelerated damage evolution, the central server generates an early warning signal for accelerated local instability of the roof. This early warning signal includes spatial positioning codes, which may include the stope number, intermediate section number, roadway number, working face number, grid cell number, three-dimensional spatial coordinate range, risk area boundary information, and control zone identifiers pre-agreed with the mine safety management system.

[0129] The central server sends early warning signals to relevant execution units via the mine safety control bus (e.g., based on CAN bus, industrial Ethernet, or intrinsically safe communication protocols for mining), triggering hardware-level interrupts or interlocking control actions. Interlocking control actions include, but are not limited to, one or more combinations of the following: stopping the operation of coal mining machines in the warning area and adjacent areas; stopping the scraper conveyor; cutting off the power supply to the designated work area; triggering audible and visual alarm devices; initiating emergency broadcasts to push personnel evacuation instructions; coordinating with the ventilation system to adjust airflow; coordinating with the drainage system to initiate emergency drainage; and uploading risk event records and current damage evolution factor distribution maps to the mine safety management platform in real time.

[0130] After all processing steps for the current monitoring period are completed, the differential reflectivity tensor calculated for the current monitoring period will be... The non-volatile storage units stored in the central server will be used for the next monitoring cycle. The differential reflectivity tensor of the previous monitoring period is obtained when executing step S3.

[0131] Therefore, this invention forms a complete feedback loop across monitoring cycles through the following data streams:

[0132] Differential reflectivity tensor of the previous monitoring period → Constructing the spatial anisotropic diffusion tensor through spectral decomposition → Reconstruct the spatial scattering sensitivity kernel function and assemble the sparse sensitivity matrix → Constrain the current period's wake interferometry inversion to obtain the relative velocity perturbation field. → Update 3D velocity model → Used for current-period reverse-time migration imaging → Current-period differential reflectivity tensor calculated via cross-period wavefield destructive solution → Store and feed back to the next cycle step S3.

[0133] This feedback closed-loop mechanism enables the present invention to continuously and dynamically track the spatiotemporal evolution of the damage structure of the roof rock mass in the mining area, and organically integrate the imaging recognition results of the previous cycle into the inversion and imaging calculation of the current cycle in the form of physical constraints, thereby realizing continuous tracking and dynamic prediction of the entire process of roof damage from initiation to expansion and accelerated instability.

[0134] Example 2

[0135] This embodiment pre-calculates weighting coefficients in the roof damage evolution factor. , The method for determining the ultimate rock mass fracture damage threshold, and the preset limit rock mass fracture damage threshold. The method for determining [the method] will be explained in further detail.

[0136] Preset weighting coefficients and The determination was made by combining indoor rock mechanics tests with in-situ calibration tests.

[0137] First, rock specimens with the same or similar lithology as the roof of the target mining area were selected for uniaxial graded loading damage evolution tests. During the test, the following data were recorded simultaneously at each loading stage: axial stress, axial strain, and circumferential strain; changes in ultrasonic longitudinal and transverse wave velocities (used to calculate the absolute value of the relative velocity disturbance). ); Reflectivity imaging results based on active ultrasonic phased array or acoustic emission (used to calculate the Frobenius norm of the differential reflectivity tensor). ), and reference damage variables calibrated using the strain equivalence assumption or the acoustic emission cumulative energy method. .

[0138] Assume the experiment yields Data samples under different damage states. Used as a reference for damage variables. To fit the target, the weighting coefficients are determined using the least squares method. and :

[0139] ;

[0140] Solving the above least squares problem yields the initial values ​​of the weighting coefficients for indoor experimental calibration. and .

[0141] Considering the differences in homogeneity and damage evolution patterns between small-scale indoor rock specimens and large-scale heterogeneous rock masses in the field, scale mapping corrections were further performed on the initial values ​​through field calibration tests in the stope. Field calibration tests can be based on one or more of the following data sources: roof fracture development areas confirmed by borehole inspection; spatial distribution records of historical roof falls or spalling events; independent acoustic detection or borehole radar detection results; or damage zoning confirmed by manual field reconnaissance. By using known damage areas in the field… and Values ​​are mapped to the corresponding damage levels. and Linear or nonlinear scaling corrections are performed to obtain the final weighting coefficients applicable to the actual mining area scale. and .

[0142] Preset ultimate rock mass fracture damage threshold The determination method is as follows: a uniaxial graded loading damage evolution test is conducted on rock specimens with the same lithology as the roof of the target monitoring area, and the roof damage evolution factor corresponding to the specimen entering the macroscopic crack penetration stage (post-peak region of the stress-strain curve) is recorded. The value is taken as the critical value and then multiplied by a preset safety reduction factor. The safety reduction factor can be selected according to the requirements of the mine safety management regulations, for example, taking the critical value. This will further improve the lead time and safety margin of early warnings.

[0143] Example 3

[0144] This embodiment provides a preferred parameter setting example to further illustrate the feasibility of implementing the present invention.

[0145] In this example, the target monitoring area is the roof of a longwall mine in an underground mine, and the monitoring cycle step is... The time was set to 30 minutes; the three-dimensional spatial mesh was divided into regular hexahedral sections with a mesh size of 2m; the initial background wave velocity field was established based on the results of on-site borehole acoustic tests.

[0146] In homology clustering, similarity threshold Set as Spatial distance threshold Set as Reference diffusion coefficient Values Minimum diffusion coefficient Pick Damage sensitivity adjustment coefficient Pick Sparse threshold of the sparse sensitivity matrix Set as .

[0147] The wake interference time window begins at three times the S-wave travel time and ends when the signal-to-noise ratio drops to twice the background noise level. Inversion regularization parameters are then retrieved. and The L-curve method is used for adaptive determination.

[0148] After calibration, the weighting coefficients Pick , Pick Ultimate rock mass fracture damage threshold Critical damage value for the macroscopic fracture penetration stage determined by indoor tests Multiply by the safety reduction factor Set later .

[0149] The interlocking control action is configured as follows: Level 1 warning ( and and ) Triggers automatic shutdown of coal mining machines and scraper conveyors in the warning area; Level II warning ( and and At the same time, the power supply to the warning area and adjacent areas is cut off, and the audible and visual alarms for the evacuation of personnel in the entire area and the emergency broadcast are activated.

[0150] The above embodiments are merely illustrative examples of the present invention and are not intended to limit the scope of protection of the present invention. Any equivalent substitutions or adaptive adjustments made by those skilled in the art to the sequence of steps, parameter selection, algorithm model, sensor deployment scheme, server hardware architecture, or early warning interlocking control strategy without departing from the technical concept of the present invention should be understood as falling within the scope of protection of the present invention.

Claims

1. A method for roof safety early warning and interlocking control based on microseismic full-waveform reverse time migration imaging, characterized in that, Includes the following steps: S1: Discretize the target monitoring area into a spatial grid, configure the initial background wave velocity field, and set the monitoring cycle sequence; S2: During the current monitoring period, the collected continuous waveform signals are dynamically truncated and clustered based on their similarity to obtain a cluster of microseismic events with similar origins. S3: Obtain the differential reflectivity tensor calculated by reverse time-shift imaging in the previous monitoring cycle, construct a spatial anisotropic diffusion tensor by performing spectral decomposition on the differential reflectivity tensor, and reconstruct the spatial scattering sensitivity kernel function based on the diffusion tensor to establish a sparse sensitivity matrix characterizing the sensitivity relationship between wave velocity disturbance and wake phase delay. S4: Using the wake sequence of the co-source microseismic cluster, combined with the sparse sensitivity matrix, inversion is performed to obtain the velocity relative disturbance field of the current monitoring period, and the three-dimensional velocity model is updated accordingly. S5: Using the updated three-dimensional velocity model as input, perform reverse time migration on the microseismic waveform of the co-source microseismic cluster, and decouple the elastic wave field during the reverse time migration process; Based on the time-difference cross-correlation imaging condition, cross-cycle wavefield cancellation is performed on the decoupled wavefield of the current monitoring period and the previous monitoring period; wherein, the cross-cycle wavefield cancellation includes matrix subtraction of the cross-correlation integral result of the current period's decoupled wavefield and the historical integral result of the previous monitoring period after event energy weighted normalization, so as to solve the differential reflectivity tensor of the current monitoring period. S6: Based on the absolute value scalar of the velocity relative to the perturbation field and the Frobenius matrix norm of the differential reflectivity tensor, a multi-source joint weighted algebraic mapping is performed to construct the roof damage evolution factor, which satisfies: ； in, The roof damage evolution factor; The relative perturbation value of wave velocity in the velocity relative perturbation field in the grid element The absolute value scalar at that point; The differential reflectivity tensor in the mesh element The Frobenius matrix norm at that location; and Preset weighting coefficients; The first-order time rate gradient and second-order time rate acceleration of the roof damage evolution factor are calculated using the following discrete difference formulas: ， ； in, The preset monitoring cycle step size; When the roof damage evolution factor and its first-order time rate gradient and second-order time rate acceleration satisfy the preset composite damage acceleration evolution criterion, an early warning signal for the accelerated evolution of local roof instability containing spatial positioning code is output, and a hardware-level interrupt or interlocking control action is triggered through the mine safety control bus. The interlocking control action includes stopping the operation of mining equipment or cutting off the power supply to the designated work area. Furthermore, the differential reflectivity tensor of the current monitoring cycle is stored as the input tensor for constructing the spatial anisotropic diffusion tensor in step S3 of the next monitoring cycle, so as to form a feedback loop across monitoring cycles.

2. The method according to claim 1, characterized in that, When the current monitoring period is the initial monitoring period, the differential reflectance tensor of the previous monitoring period is initialized to a preset three-dimensional all-zero tensor, or initialized to a reference reflectance tensor obtained based on the complete roof rock mass background model.

3. The method according to claim 1, characterized in that, In step S3, a spatially anisotropic diffusion tensor is constructed by performing spectral decomposition on the differential reflectivity tensor, specifically as follows: The differential reflectivity tensor, which is a three-dimensional second-order symmetric tensor, is subjected to eigenvalue decomposition to obtain eigenvalue scalars and corresponding unit eigenvectors; Construct the corresponding directional projection tensor based on the outer product of each unit eigenvector and its transpose, and perform a weighted summation of each directional projection tensor according to the corresponding eigenvalues. Then, perform an algebraic summation of the summation result with the three-dimensional unit matrix scaled by the preset background isotropic diffusion constant to obtain the spatial anisotropic diffusion tensor.

4. The method according to claim 1, characterized in that, The event energy weighted normalization process includes: If there are no historical events corresponding to the current microseismic cluster in the previous monitoring period, then the all-zero tensor is used as the historical integration result. If multiple historical events corresponding to the current co-source microseismic cluster exist within the previous monitoring period, then the historical conjugate wavefield tensor integral results of each event are weighted and averaged using the scalar seismic moments of each event as weights to obtain the standardized historical integral results.

5. The method according to claim 1, characterized in that, The preset weighting coefficient and The data were obtained through uniaxial graded loading damage evolution tests on homogeneous rock specimens combined with least squares fitting calibration, and then corrected by scale mapping through on-site calibration tests in the mining area.

6. The method according to claim 1, characterized in that, The criterion for accelerated evolution of composite damage is: ； in, To preset the ultimate rock mass fracture damage threshold, the The critical damage value for macroscopic fracture penetration was determined based on uniaxial graded loading damage evolution tests on homogeneous rock specimens.

7. The method according to claim 1, characterized in that, The updating of the three-dimensional velocity model in step S4 specifically involves: recursively updating the three-dimensional velocity model of the previous monitoring period based on the velocity relative disturbance field, so that the three-dimensional velocity model evolves continuously in the monitoring period sequence.

8. The method according to claim 1, characterized in that, The inversion in step S4 is solved using regularization constraints, which include at least one of spatial smoothing constraints, Tikhonov regularization constraints, or anisotropy constraints based on geological structure orientation.

Images