A biomimetic rock-breaking tool load adaptive control method and system

Abstract

This invention discloses a biomimetic rock-breaking tool load adaptive control method and system, belonging to the field of load adaptive control technology. The method specifically includes: collecting a dataset of the biomimetic rock-breaking tool's operation during continuous rock-breaking; based on the dataset, performing phase organization of the high-load process to obtain a high-load rock-breaking process sequence; based on the high-load rock-breaking process sequence, constructing a high-load process phase correlation diagram, and performing recursive estimation of the switching state space under the closed-loop lag constraint of the rock-breaking result to obtain a high-load state result set; based on the high-load state result set, performing rolling prediction of the next rock-breaking phase state transition and closed-loop correction of the control quantity to obtain a load adaptive control instruction set. This application improves the accuracy and stability of the biomimetic rock-breaking tool load adaptive control.

Description

Technical Field

[0001] This invention belongs to the field of load adaptive control, specifically relating to a biomimetic rock-breaking tool load adaptive control method and system. Background Technology

[0002] In the continuous rock-breaking process of biomimetic rock-breaking tools, load adaptive control needs to determine the current rock-breaking state based on the tool load and rock-breaking response, and adjust the tool speed and feed parameters accordingly. Biomimetic rock-breaking tools typically apply local forces to the rock mass using ridges, curved surfaces, or non-smooth structures. However, crack initiation, crack propagation, and tool obstruction can all cause high loads, resulting in mixed load characteristics under different rock-breaking states, making it difficult to accurately set control parameters.

[0003] In existing technologies, biomimetic rock-breaking tool load adaptive control methods typically adjust the tool rotation speed and feed parameters based on changes in the total tool torque, axial load, and drive current. However, the short-term high load generated during effective fracturing and the high load generated by ineffective tool obstruction exhibit similar behavior. Existing methods do not consider load release, vibration response, and phase changes in adjacent rock-breaking cycles after high load to distinguish between the high load generated during effective fracturing and the high load generated during ineffective obstruction. This leads to the effective fracturing process being easily misjudged as an ineffective obstruction process, resulting in an incorrect reduction in tool rotation speed or feed rate. Consequently, the accuracy and stability of biomimetic rock-breaking tool load adaptive control are poor. Summary of the Invention The purpose of this invention is to provide a biomimetic rock-breaking tool load adaptive control method and system. By constructing a phase correlation diagram of the high-load process and performing recursive estimation of the switching state space under the closed lag constraint of the rock-breaking result, the effective fracturing high load and the ineffective obstructed high load are distinguished, thereby improving the accuracy and stability of the biomimetic rock-breaking tool load adaptive control.

[0004] To solve the above-mentioned technical problems, the present invention provides the following technical solution: On one hand, the present invention provides a biomimetic rock-breaking tool load adaptive control method, comprising: The operation dataset of the biomimetic rock-breaking tool during the continuous rock-breaking process is collected. The operation dataset includes load data, motion position data, tool rotation speed, feed speed, and control execution data. Based on the running dataset, phase organization of the high-load process is performed to obtain the high-load rock breaking process sequence; Based on the high-load rock breaking process sequence, a phase correlation diagram of the high-load process is constructed, and a recursive estimation of the switching state space under the closed lag constraint of the rock breaking result is performed to obtain a high-load state result set, which includes an effective fracturing state sequence, an ineffective obstructed state sequence, and a high-load state transition response relationship. Based on the high-load state result set, the rolling prediction of the next rock-breaking phase state transition and the closed-loop correction of the control quantity are performed to obtain the load adaptive control instruction set.

[0005] Specifically, based on the runtime dataset, high-load process phase organization is performed to obtain a high-load rock-breaking process sequence, including: Based on motion position data, the continuous rock-breaking process is divided into rock-breaking cycles, and the rotation phase within each rock-breaking cycle is determined to obtain the rock-breaking phase sequence. Based on the rock-breaking phase sequence, the in-phase load benchmark of the load data under each rotation phase is determined, and the deviation sequence is calculated; Based on the deviation sequence, high-load sections are extracted, and the rock-breaking cycle and rotation phase range corresponding to each high-load section are determined. Based on the rock-breaking cycle and rotation phase range corresponding to each high-load section, in-phase data of adjacent rock-breaking cycles are extracted and correlated with each high-load section to obtain the high-load rock-breaking process sequence.

[0006] Specifically, based on the high-load rock-breaking process sequence, a phase correlation diagram of the high-load process is constructed, and a recursive estimation of the switching state space under the closed lag constraint of the rock-breaking results is performed to obtain the high-load state result set, including: Based on the high-load rock breaking process sequence, the deviation sequence change corresponding to each high-load section was determined, and the load rising process, peak action process and post-peak response process were divided to obtain the process stage division results. Based on the high-load rock breaking process sequence and process stage division results, phase synchronization mode extraction was performed, and process stage correlation was performed to obtain high-load mode fragments; Based on high-load modal fragments, cross-channel homology aggregation is performed to obtain a set of high-load homology fragments; Construct a phase correlation graph for high-load processes based on a set of high-load homogeneous fragments; Based on the phase correlation diagram of the high-load process, a recursive estimation of the switching state space under the closed lag constraint of the rock breaking result is performed to obtain the high-load state result set.

[0007] Specifically, based on high-load modal fragments, cross-channel homology aggregation is performed to obtain a set of high-load homology fragments, including: Based on the high-load modal segments, the central order and modal phase are extracted, the signal channel identifier corresponding to each high-load modal segment is determined, and the rotation phase range is associated to obtain the modal segment feature set; Based on the modal segment feature set, high-load modal segments with the same center order, overlapping rotation phase range, and different signal channel identifiers are combined to obtain cross-channel candidate segment groups; Based on the cross-channel candidate fragment group, the response phase difference of each cross-channel candidate fragment within the overlapping rotation phase range is calculated, and the degree of fluctuation of the response phase difference is determined. Based on the degree of fluctuation of the response phase difference, cross-stage stable relationships are determined, and homogeneous fragment aggregation and attribute merging are performed to obtain a high-load homogeneous fragment set.

[0008] Specifically, based on a set of high-load homogeneous fragments, a phase correlation graph of the high-load process is constructed, including: Based on the set of high-load homologous fragments, each high-load homologous fragment is mapped to a high-load action node, and associated with the corresponding high-load section, rock breaking cycle, rotation phase range and process stage to obtain the set of high-load action nodes. Based on the set of high-load action nodes, the node stage succession relationship within the same high-load segment is determined, process evolution edges are established, and high-load formation nodes, peak action nodes, and post-peak response nodes are identified. Based on the rock-breaking cycle and rotation phase range of the post-peak response node, corresponding data are matched from the in-phase data of adjacent rock-breaking cycles, mapped to in-phase verification nodes, and result verification edges are established. Based on the control execution data, determine the execution response range and the corresponding high-load action node, and establish the control action edge; Based on the high-load action node set, in-phase verification nodes, process evolution edges, result verification edges, and control action edges, a graph structure is organized to obtain a high-load process phase correlation graph.

[0009] Specifically, based on the phase correlation diagram of the high-load process, a recursive estimation of the switching state space under the closed lag constraint of the rock breaking results is performed to obtain the high-load state result set, including: Based on the phase correlation graph of the high-load process, the result verification path connecting the high-load forming node and the same-phase verification node is extracted; The high-load forming node is taken as the state node to be confirmed, and a switching state space model is constructed by combining the high-load homogeneous fragment set and control execution data. Based on the result verification path and the switching state space model, a forward state recursive estimation is performed on the state node to be confirmed to obtain the forward state estimation result. Based on the forward state estimation results, closed evidence-driven backward state smoothing and control response correlation are performed to obtain a high-load state result set.

[0010] Specifically, based on the result verification path and the switching state space model, a forward recursive estimation of the state node to be confirmed is performed to obtain the forward state estimation result, including: Based on the process evolution edge, the high load forming nodes, peak effect nodes, and post-peak response nodes are sequentially organized to obtain the node evolution sequence; Based on the node evolution sequence and control action edge, the control execution data is associated with the corresponding high-load action node to obtain the node observation input sequence; Based on the switching state space model, node evolution sequence and node state input sequence, forward state recursive estimation is performed to obtain node state probability sequence and state transition probability sequence. Based on the node state probability sequence and the state transition probability sequence, the node is kept in a state of unconfirmation under high load before the result verification path is closed, thus obtaining the forward state estimation result.

[0011] Specifically, based on the forward state estimation results, closed-loop evidence-driven backward state smoothing and control response correlation are performed to obtain a high-load state result set, including: Based on the forward state estimation results and result verification path, the state smoothing interval is determined and the probability sequence is organized to obtain the state smoothing interval and the probability sequence to be smoothed. Based on the state smoothing interval, the post-peak process response change and the in-phase load change of adjacent periods are calculated to obtain the closed-loop verification feature sequence. Based on the closure verification feature sequence, valid fracturing closure evidence and invalid blocked closure evidence are identified, and the closure evidence results are obtained. Based on the probability sequence to be smoothed and the closure evidence results, reverse state smoothing is performed along the result verification path to obtain the closed state probability sequence and the closed state transition probability sequence. Based on the closed-state probability sequence and the closed-state transition probability sequence, high-load state determination and control response correlation are performed to obtain a high-load state result set.

[0012] Specifically, based on the high-load state result set, the next rock-breaking phase state transition rolling prediction and control quantity closed-loop correction are performed to obtain the load adaptive control instruction set, including: Based on the effective fracturing state sequence, the ineffective obstructed state sequence, the tool rotation speed, the feed speed and the control execution data, the current high load state and the current control state are determined, and the current state information is obtained. Based on the current state information, a set of candidate control combinations for tool rotation speed and feed rate is generated; Based on the current high load state, candidate control combination set and high load state transition response relationship, the state transition rolling prediction of the next rock breaking phase is carried out to obtain the candidate state prediction result set; Based on the current control state and the candidate state prediction result set, the target control combination is determined, and the control quantity is closed-loop corrected to obtain the load adaptive control instruction set.

[0013] On the other hand, the present invention provides a biomimetic rock-breaking tool load adaptive control system, comprising: The data acquisition module is used to collect the operational dataset of the bionic rock-breaking tool during the continuous rock-breaking process; The phase organization module performs phase organization of the high-load process based on the running dataset to obtain the high-load rock breaking process sequence; The recursive estimation module constructs a phase correlation diagram of the high-load process based on the high-load rock breaking process sequence, and performs recursive estimation of the switching state space under the closed lag constraint of the rock breaking result to obtain the high-load state result set. The prediction and correction module, based on the high load state result set, performs rolling prediction of the next rock breaking phase state transition and closed-loop correction of the control quantity to obtain the load adaptive control instruction set.

[0014] Compared with the prior art, the beneficial effects of the present invention include: by constructing a phase correlation diagram of the high-load process that includes process evolution relationship, result verification relationship and control action relationship, the high-load formation process is correlated with the post-peak response and the in-phase change of adjacent rock breaking cycles, thereby improving the accuracy of distinguishing between effective fracturing high load and ineffective obstructed high load; by using the switching state space recursive estimation under the closed hysteresis constraint of rock breaking result, the forward state estimation result is corrected using the closed evidence after the formation of high load, and the tool speed and feed parameters are adjusted differentially, thereby enhancing the stability of load adaptive control. Attached Figure Description

[0015] Figure 1 A flowchart of a biomimetic rock-breaking tool load adaptive control method provided by the present invention; Figure 2 This is a schematic diagram of the closure verification feature acquisition provided by the present invention; Figure 3 This is a schematic diagram of the closed evidence determination and fusion provided by the present invention; Figure 4 This is a schematic diagram illustrating the correlation between reverse state smoothing and control provided by the present invention; Figure 5 The present invention provides a structural diagram of a biomimetic rock-breaking tool load adaptive control system. Detailed Implementation

[0016] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the embodiments of the present invention and the specific features in the embodiments are detailed descriptions of the technical solution of the present invention, rather than limitations thereof. In the absence of conflict, the embodiments of the present invention and the technical features in the embodiments can be combined with each other.

[0017] Example 1: See Figures 1-4 This embodiment provides a biomimetic rock-breaking tool load adaptive control method, including the following specific steps: Step S1: Collect the operational dataset of the biomimetic rock-breaking tool during the continuous rock-breaking process.

[0018] In this embodiment, the continuous rock breaking process is defined as the time interval during which both the rotary actuator and the propulsion actuator are in operation and the biomimetic rock-breaking tool is in continuous contact with the rock mass.

[0019] The total torque of the tool is collected by a torque sensor installed on the tool spindle drive path, the axial load is collected by an axial load sensor installed on the tool feed force path, and the drive current is collected by a current detection unit of the rotary drive device and the feed drive device to obtain the raw load data.

[0020] Static sampling data is acquired when the cutting tool is not in contact with the rock mass and the rotary and propulsion actuators are not applying rock-breaking action. Based on this static sampling data, the zero-point offset of each load channel is determined. Specifically, the method is as follows: all static sampling values ​​for the same load channel are summed, and the summation result is divided by the number of static sampling values ​​to obtain the zero-point offset of the corresponding load channel. The corresponding zero-point offset is then subtracted from the load sampling values ​​during continuous rock-breaking to obtain the zero-point correction result.

[0021] The dimension conversion factor is determined based on the change in standard load and the change in sensor output applied during sensor calibration. Specifically, the method is as follows: divide the change in standard load by the corresponding change in sensor output to obtain the dimension conversion factor. Multiply the zero-point correction result by the corresponding dimension conversion factor to obtain the load data.

[0022] The tool rotation angle is acquired using a rotary encoder mounted on the tool spindle, and the propulsion displacement is acquired using a displacement sensor mounted on the propulsion mechanism. The specified rotation direction of the tool is defined as the positive direction of the rotation angle, and the direction in which the tool enters the rock mass is defined as the positive direction of the propulsion displacement. By unifying the directions of the rotation angle and the propulsion displacement, the motion position data is obtained.

[0023] Based on the speed feedback terminals of the rotary actuator and the feed actuator, the tool rotation speed and feed speed are collected respectively. The tool rotation speed direction is verified according to the direction of rotation angle change, and the feed speed direction is verified according to the direction of feed displacement change. Sampling intervals where the speed direction and the corresponding position change direction are inconsistent are marked as speed abnormal intervals.

[0024] The controller, based on a biomimetic rock-breaking tool, collects rotation control commands, propulsion control commands, and feedback data from the actuators. The rotation control commands include the target tool speed and rotational operating status; the propulsion control commands include the target propulsion speed and propulsion operating status; and the actuator feedback data includes the actual tool speed, actual propulsion speed, rotation tracking deviation, propulsion tracking deviation, rotation execution status, and propulsion execution status. Control execution data is obtained by performing time-series correlation based on the controller output time.

[0025] A unified sampling time sequence is established based on the controller system clock, and load data, motion position data, tool speed, feed speed and control execution data are associated with the corresponding sampling time.

[0026] The allowable completion time is set based on the interpolation error of complete running data under simulated missing conditions. The specific setting method is as follows: Continuous sampling points are sequentially deleted from historical complete running data, and linear interpolation is performed based on the effective sampling values ​​before and after the missing interval. The original values ​​of the deleted sampling points are subtracted from the interpolation result, and the absolute value of the difference is taken to obtain the interpolation error. The interpolation error corresponding to each consecutive missing length is compared with three times the corresponding sensor repeatability error. The maximum number of consecutive missing sampling points when the interpolation error does not exceed three times the corresponding sensor repeatability error is multiplied by the sampling period to obtain the allowable completion time.

[0027] Linear interpolation is performed on data with consecutive missing durations not exceeding the allowed completion time, and data intervals with consecutive missing durations exceeding the allowed completion time are marked as invalid data intervals.

[0028] Based on a unified sampling time sequence, the load data, motion position data, tool speed, feed speed, control execution data, speed abnormal intervals, and invalid data intervals are arranged in time sequence to obtain the running dataset.

[0029] Step S2: Based on the running dataset, perform high-load process phase organization to obtain the high-load rock breaking process sequence.

[0030] The specific steps of step S2 are as follows: Step S201: Based on the motion position data, divide the continuous rock breaking process into rock breaking cycles and determine the rotation phase within each rock breaking cycle to obtain the rock breaking phase sequence.

[0031] In this embodiment, the number of repeating biomimetic rock-breaking structures along the circumference of the tool is read based on the tool's structural parameters. These repeating biomimetic rock-breaking structures are ridges, protrusions, curved surface units, or non-smooth structures that periodically and repeatedly act on the rock mass during tool rotation. The basic rock-breaking cycle angle is obtained by dividing the rotation angle corresponding to one revolution by the number of repeating biomimetic rock-breaking structures. When the tool does not have a circumferential repeating structure, the rotation angle of one revolution is determined as the basic rock-breaking cycle angle.

[0032] Angle wrapping is determined based on the tool rotation angle at adjacent sampling times. When the subsequent rotation angle is less than the previous rotation angle, and the difference between the previous and subsequent rotation angles exceeds half a rotation angle, the cumulative rotation count is increased by one revolution. The half-rotation angle is obtained by dividing the full rotation angle by 2 and is used to distinguish between normal reverse fluctuations and angle wrapping. The unfolded rotation angle is obtained by adding the current rotation angle and the integer rotation angle corresponding to the cumulative full rotation count.

[0033] Subtract the rotation angle at the start of the continuous rock-breaking process from the rotation angle at each sampling time to obtain the cumulative change in rotation angle. Divide the cumulative change in rotation angle by the basic rock-breaking cycle angle and round down to obtain the rock-breaking cycle number. Multiply the rock-breaking cycle number by the basic rock-breaking cycle angle to obtain the cumulative cycle angle corresponding to the completed rock-breaking cycle. Subtract the cumulative cycle angle from the cumulative change in rotation angle to obtain the rotation phase angle within the cycle.

[0034] The rotation phase division step size is set based on the minimum resolution angle of the rotary encoder and the actual rotation sampling interval during continuous rock breaking. The specific setting method is as follows: subtract the unfolded rotation angle of the previous valid sampling moment from the unfolded rotation angle of the subsequent valid sampling moment to obtain the angle difference between adjacent sampling moments. Then, arrange all adjacent angle differences in ascending order and extract the value corresponding to the middle position to obtain the median of the adjacent angle differences. The larger value between this median and the minimum resolution angle of the rotary encoder is taken as the rotation phase division step size.

[0035] Divide the basic rock-breaking cycle angle by the rotation phase step size and round down to obtain the number of rotation phases within each rock-breaking cycle. Divide the rotation phase angle within the cycle by the rotation phase step size and round down to obtain the rotation phase number. Correlate the rock-breaking cycle number, rotation phase number, rotation phase angle, sampling time, and advance displacement to obtain the rock-breaking phase sequence.

[0036] Step S202: Based on the rock-breaking phase sequence, determine the in-phase load reference of the load data under each rotation phase, and calculate the deviation sequence.

[0037] In this embodiment, based on the rock-breaking phase sequence, the total tool torque, axial load, and drive current are associated with the corresponding rock-breaking cycle number and rotation phase number to obtain phased load data.

[0038] The in-phase reference window consists of consecutive and valid complete rock-breaking cycles preceding the current rock-breaking cycle. The length of the in-phase reference window is set based on the stability of the in-phase load statistics. Specifically, the setting method is as follows: starting from the current rock-breaking cycle, increment by one complete rock-breaking cycle. For each additional complete rock-breaking cycle, recalculate the median load for each rotational phase. Subtract the median load before the extended window from the median load after the extended window, and take the absolute value of the difference to obtain the load reference change.

[0039] The load reference stability threshold is set based on the repeatability error of the corresponding sensor, specifically three times the repeatability error. The repeatability error is determined by measuring values ​​obtained from applying the same standard load multiple times during sensor calibration. The specific method for determining the repeatability error is as follows: sum all measured values, divide the sum by the number of measured values ​​to obtain the average value; subtract the average value from each measured value to obtain the measurement deviation for each value; square each measurement deviation, sum the squares of all measurement deviations, and divide by the number of measured values ​​minus 1 to obtain the mean square of the measurement deviations; take the square root of the mean square of the measurement deviations to obtain the repeatability error.

[0040] The number of consecutive stability determinations is set based on the longest consecutive period in which the load baseline change exceeds the load baseline stability threshold in historical stable rock breaking data. The specific method is as follows: count the maximum number of consecutive periods in which the load baseline change exceeds the load baseline stability threshold in historical stable rock breaking data, and then increment this maximum number of periods by 1 to obtain the number of consecutive stability determinations.

[0041] When the load reference change of each rotation phase continuously reaches the number of consecutive stability determinations and does not exceed the load reference stability threshold, the addition of complete rock breaking cycles is stopped, and the number of complete rock breaking cycles contained in the current window is determined as the length of the same phase reference window.

[0042] Based on the same-phase reference window, the total tool torque, axial load, and drive current are collected according to the rotational phase sequence. The load data of the same type under the same rotational phase are arranged in ascending order, and the load value corresponding to the middle position is extracted to obtain the corresponding same-phase load reference.

[0043] The discrete scale of the in-phase load is determined based on the in-phase load reference. The specific calculation process is as follows: subtract the corresponding in-phase load reference from each load sample value, take the absolute value of the difference, arrange all absolute differences in ascending order, and extract the absolute difference corresponding to the middle position to obtain the discrete scale of the in-phase load. When the discrete scale of the in-phase load is less than the minimum resolution value of the corresponding sensor, the minimum resolution value of the sensor is determined as the discrete scale of the in-phase load.

[0044] The positive load difference is obtained by subtracting the corresponding in-phase load reference from the current load sample value. When the positive load difference is less than zero, the standardized positive deviation is set to zero. When the positive load difference is not less than zero, the positive load difference is divided by the corresponding in-phase load discrete scale to obtain the standardized positive deviation.

[0045] The combined deviation is calculated based on the standardized positive deviations of the total tool torque, axial load, and drive current. The specific calculation process is as follows: square the three standardized positive deviations separately, add the three squared results, divide the sum by the number of load channels, and take the square root of the result to obtain the combined deviation. The standardized positive deviations and combined deviations of each channel are then organized according to a unified sampling time to obtain the deviation sequence.

[0046] Step S203: Based on the deviation sequence, extract the high-load segments and determine the rock breaking cycle and rotation phase range corresponding to each high-load segment.

[0047] In this embodiment, the high load entry threshold is set based on the allowable false triggering ratio of the combined deviation under normal rock breaking conditions. Specifically, the setting method is as follows: divide the allowable false triggering sampling points by the total number of normal rock breaking sampling points to obtain the false triggering ratio. Arrange the combined deviations within the in-phase reference window in ascending order, extract the combined deviation corresponding to the point where the cumulative ratio reaches one, and subtract the false triggering ratio to obtain the high load entry threshold.

[0048] The high-load exit threshold is set based on the allowable repeated switching ratio during the high-load exit process. Specifically, the method is as follows: candidate exit thresholds are selected sequentially from candidate joint deviations that are less than the high-load entry threshold. Exit decisions are made for historical high-load segments based on each candidate exit threshold, and the number of repeated switching cycles between entry and exit during the high-load state is counted. The candidate exit threshold with the largest number of repeated switching cycles not exceeding the allowable number of repeated switching cycles is determined as the high-load exit threshold, ensuring that the high-load exit threshold is less than the high-load entry threshold.

[0049] The auxiliary judgment thresholds for each load channel are set according to the allowable false triggering ratio of the standardized positive deviation under normal rock breaking conditions. The specific setting method is the same as that for the high load entry threshold, resulting in the auxiliary judgment thresholds for total tool torque, axial load, and drive current.

[0050] The minimum number of continuous sampling points under high load is set based on the effective bandwidth of the load sensor and the duration of transient interference. The specific setting method is as follows: Take the reciprocal of the lowest effective bandwidth among the tool total torque sensor, axial load sensor, and drive current detection unit to obtain the shortest reliable response time. Divide the shortest reliable response time by the load data sampling period and round the result up to obtain the bandwidth-constrained sampling point number. Count the longest continuous sampling point for single-point pulse interference during sensor calibration, and then increase this longest continuous sampling point number by 1 to obtain the anti-interference sampling point number. Take the larger value between the bandwidth-constrained sampling point number and the anti-interference sampling point number to obtain the final minimum number of continuous sampling points.

[0051] The sampling time at which the combined deviation is not less than the high load entry threshold, and the standardized positive deviation of at least two load channels is not less than the corresponding channel's auxiliary judgment threshold, is marked as a high load entry candidate sampling time. When the number of consecutive high load entry candidate sampling times reaches the final minimum number of consecutive sampling points, the first candidate sampling time is determined as the start time of the high load segment.

[0052] Traverse the deviation sequence backward from the start time of the high-load segment. When the number of sampling points where the joint deviation is continuously lower than the high-load exit threshold reaches the final minimum number of consecutive sampling points, the last sampling time before that consecutive sampling interval is determined as the end time of the high-load segment.

[0053] The maximum mergeable interval for high-load segments is set based on the shortest reliable response time. Specifically, the shortest reliable response time is determined as the maximum mergeable interval. When the time interval between two adjacent high-load segments is not greater than the maximum mergeable interval, the period from the start time of the preceding high-load segment to the end time of the following high-load segment is determined as the merged high-load segment.

[0054] Based on the rock-breaking phase sequence, the rock-breaking cycle number and rotation phase number corresponding to the start and end times of each high-load segment are read to obtain the rock-breaking cycle and rotation phase range corresponding to each high-load segment. For high-load segments that cross the rock-breaking cycle boundary, the phase range at the end of the cycle and the phase range at the beginning of the next cycle are recorded respectively, and they are associated as the same high-load segment.

[0055] Step S204: Based on the rock breaking cycle and rotation phase range corresponding to each high-load section, extract the in-phase data of adjacent rock breaking cycles and associate them with each high-load section to obtain the high-load rock breaking process sequence.

[0056] In this embodiment, the preceding complete rock-breaking cycle of the current high-load section is defined as the preceding adjacent rock-breaking cycle, and the following adjacent complete rock-breaking cycle is defined as the following adjacent rock-breaking cycle.

[0057] Based on the starting and ending rotation phases of the high-load section, load data, motion position data, tool speed, feed speed, and control execution data that fall within the same rotation phase range are extracted from the preceding adjacent rock breaking cycle and the following adjacent rock breaking cycle, respectively, to obtain the same-phase data of the preceding cycle and the same-phase data of the following cycle.

[0058] When the high-load section crosses the boundary of the rock breaking cycle, the same-phase data at the end of the cycle and the beginning of the next cycle are extracted respectively, and then spliced ​​according to the rotation phase order.

[0059] If the sampling positions of the current high-load segment and the adjacent period's same-phase data are inconsistent, a unified rotating phase position is established according to the rotating phase division step size determined in step S201. For any target rotating phase position, its two preceding and following effective sampling positions are determined. The forward phase distance is obtained by subtracting the preceding effective sampling position from the target rotating phase position. The backward phase distance is obtained by subtracting the target rotating phase position from the following effective sampling position. The interpolation weight of the preceding effective sample value is obtained by dividing the backward phase distance by the sum of the preceding and backward phase distances. The interpolation weight of the following effective sample value is obtained by dividing the forward phase distance by the sum of the two phase distances. The two effective sample values ​​are multiplied by their corresponding interpolation weights, and then the two products are added together to obtain the data value of the target rotating phase position.

[0060] The event sequence is obtained by assigning event numbers according to the occurrence time of the high-load section, associating the event number, rock breaking cycle, rotation phase range, current high-load section data, previous cycle same-phase data, subsequent cycle same-phase data and corresponding control execution data.

[0061] Step S3: Based on the high-load rock breaking process sequence, construct the phase correlation diagram of the high-load process, and perform recursive estimation of the switching state space under the closed lag constraint of the rock breaking result to obtain the high-load state result set.

[0062] The specific steps of step S3 are as follows: Step S301: Based on the high-load rock breaking process sequence, determine the deviation sequence change corresponding to each high-load section, divide the load rising process, peak action process and post-peak response process, and obtain the process stage division results.

[0063] In this embodiment, the joint deviation sequence of each high-load segment is extracted according to the event sequence number, and the joint deviation sequence is smoothed by median filtering.

[0064] The median filter window width is set according to the final minimum number of consecutive sampling points determined in step S203. Specifically, the final minimum number of consecutive sampling points is used as the initial window width. When the initial window width is even, it is increased by 1 so that the filter window contains an odd number of sampling points, thus obtaining the median filter window width.

[0065] Extract the sampling time corresponding to the maximum value of the smoothed deviation sequence to obtain the peak sampling time. The peak effect determination line is determined based on the upper quartile statistical position of the smoothed deviation in the high-load segment. The specific determination method is as follows: arrange the smoothed deviations in the high-load segment in ascending order, multiply the number of smoothed deviations by three-quarters to obtain the upper quartile statistical position, extract the smoothed deviation corresponding to that position, and obtain the peak effect determination line.

[0066] The peak action process is obtained by finding consecutive sampling points that are not lower than the peak action determination line both forward and backward from the peak sampling time. The interval from the start of the high load section to the beginning of the peak action process is defined as the load rise process, and the interval from the end of the peak action process to the end of the high load section is defined as the post-peak response process.

[0067] When the number of sampling points in the load rise process or the post-peak response process is less than the final minimum number of continuous sampling points determined in step S203, the corresponding process is incorporated into the peak action process. Each process is then correlated with its corresponding rotational phase range to obtain the process stage division results.

[0068] Step S302: Based on the high-load rock breaking process sequence and process stage division results, phase synchronization mode extraction is performed, and process stage correlation is performed to obtain high-load mode fragments.

[0069] In this embodiment, based on the high-load rock breaking process sequence, the total tool torque, axial load, and drive current in each high-load section are extracted according to the high-load section identifier.

[0070] Based on the rotation phase division step size determined in step S201, equally spaced rotation phase positions are established within the rotation phase range corresponding to each high-load section.

[0071] Using the rotational phase interpolation method in step S204, the total tool torque, axial load, and drive current are mapped to each equally spaced rotational phase position to obtain the phase sequence of the total tool torque, the phase sequence of the axial load, and the phase sequence of the drive current.

[0072] Based on the in-phase load reference and in-phase load discrete scale obtained in step S202, the phase sequence of each load channel is standardized. For any load channel and any rotational phase position, the corresponding in-phase load reference is subtracted from the load value corresponding to that rotational phase position to obtain the in-phase load difference. The in-phase load difference is divided by the corresponding in-phase load discrete scale to obtain the standardized phase value of the load channel at that rotational phase position. The standardized phase values ​​are arranged according to the rotational phase position to obtain a multi-channel standardized phase sequence, including the standardized phase sequence of total tool torque, the standardized phase sequence of axial load, and the standardized phase sequence of drive current.

[0073] Multivariate variational mode decomposition (MMD) is employed to extract phase-synchronous modes from multi-channel standardized phase sequences. Conventional MMD typically organizes channel signals according to sampling time and shares the center frequency for the same mode number across different channels. However, during continuous rock breaking by a biomimetic rock-breaking tool, the tool rotation speed changes, and the frequency of the signal corresponding to the same rock-breaking action on the time axis varies with the tool rotation speed. If decomposition is performed directly according to time and frequency, the total torque response, axial load response, and drive current response belonging to the same rotational order may be classified into different modes. To address this issue, this application adjusts the independent variable of the MMD from sampling time to rotational phase and adjusts the center frequency of the same mode number in different load channels to a shared center order. The shared center order represents the number of times the corresponding mode repeats within a complete rotational cycle, ensuring that load responses with the same rotational order are still organized into the same mode number when the tool rotation speed changes.

[0074] The constraint objective of multivariable variational mode decomposition is to minimize the sum of the phase change bandwidths of each mode after removing the corresponding shared center order, while ensuring that the superposition result of all modes in the same load channel can reconstruct the normalized phase sequence of that load channel. Specifically, for any candidate mode, a Hilbert transform is first performed on the candidate mode to obtain an analytical mode composed of the original candidate mode and its orthogonal components. Based on the current shared center order, the analytical mode is shifted to move its center order to near the zero order. The change in analytical mode after the order shift is calculated along the rotation phase direction, and the change at each rotation phase position is squared. The squared results of the same candidate mode at all rotation phase positions are summed to obtain the bandwidth of the candidate mode. In the same way, the bandwidths of all candidate modes in the tool total torque channel, axial load channel, and drive current channel are obtained. All bandwidths are summed to obtain the total bandwidth of the multi-channel modes.

[0075] Simultaneously, all candidate modes of the same load channel are superimposed according to their rotated phase positions. The normalized phase sequence of the load channel is subtracted from the superposition result to obtain the signal reconstruction residual of the load channel. The signal reconstruction residuals are squared and summed to obtain the signal reconstruction error. The optimization objective of multivariable variational mode decomposition is jointly established by the constraint of minimizing the sum of the multi-channel mode bandwidth and the constraint of the signal reconstruction error approaching zero.

[0076] The optimization objective of multivariable variational mode decomposition is solved iteratively using the alternating direction multiplier method. The parameters to be solved by the alternating direction multiplier method include the candidate modes of each load channel, the order of the shared center corresponding to each mode number, and the reconstruction constraint multipliers corresponding to each load channel.

[0077] Before the first iteration, each candidate mode is initialized to a zero-value phase sequence. The initial value of the shared center order is set at equal intervals based on the order range that can be resolved in the current high-load segment. Specifically, the zero order to the maximum resolvable order under rotating phase sampling conditions is divided into intervals equal to the number of candidate modes, and the middle order of each interval is sequentially determined as the initial shared center order of each candidate mode. The reconstruction constraint multipliers corresponding to each load channel are initialized to a zero-value sequence.

[0078] In any iteration, candidate modes are updated sequentially according to their modality numbers. When updating the current candidate mode, the phase sequences of other candidate modes that have already been updated in this round are superimposed according to the same rotational phase position to obtain the superimposed sequence of updated modes in this round. The superimposed sequence of updated modes in this round is subtracted from the normalized phase sequence of the current load channel to obtain the first reconstruction remainder sequence. Then, the phase sequences of other candidate modes that have not yet been updated in this round are superimposed according to the same rotational phase position in the previous iteration to obtain the superimposed sequence of modes that have not been updated in the previous round. The superimposed sequence of modes that have not been updated in the previous round is subtracted from the first reconstruction remainder sequence to obtain the reconstruction remainder corresponding to the current candidate mode.

[0079] Divide the reconstruction constraint multiplier of the previous round of the current load channel by 2, and then add it to the reconstruction residue corresponding to the current candidate mode to obtain the residue to be decomposed for the current candidate mode. Perform a Discrete Fourier Transform on the residue to be decomposed to obtain the frequency domain value of the residue at each discrete order position.

[0080] For any discrete order position, calculate the order difference between that discrete order position and the shared center order of the current candidate mode. Square the order difference, multiply it by 2 and the bandwidth penalty parameter to obtain the order deviation penalty. Add 1 to the order deviation penalty to obtain the mode update denominator for that discrete order position. Divide the frequency domain value of the remaining quantity to be decomposed at that discrete order position by the mode update denominator to obtain the update value of the current candidate mode at that discrete order position. Organize the update values ​​according to all discrete order positions to obtain the update frequency domain sequence of the current candidate mode.

[0081] In the same manner, the frequency domain sequences of the current candidate modes in the tool total torque channel, axial load channel, and drive current channel are updated respectively.

[0082] After updating the same mode number across all load channels, update the shared center order corresponding to that mode number. When updating the shared center order, take the amplitude of the frequency domain value of the mode at any load channel and any positive order position, and square the amplitude to obtain the mode energy at that positive order position in that load channel. Multiply the order value at that positive order position by the corresponding mode energy to obtain the order energy product at that positive order position. Add the order energy products of the same mode across all load channels and all positive order positions to obtain the weighted sum of order energies. Add the mode energies of the same mode across all load channels and all positive order positions to obtain the total modal energy. Divide the weighted sum of order energies by the total modal energy to obtain the updated shared center order.

[0083] When the sum of modal energies is less than the lower limit of modal energy, the order of the shared center of the previous round of that mode remains unchanged. The lower limit of modal energy is set according to the minimum non-zero energy in the normalized phase sequence, specifically: the minimum non-zero square value of each normalized phase value is determined as the lower limit of modal energy.

[0084] After updating all candidate modes and the shared center order, the reconstruction constraint multipliers for each load channel are updated. When updating the reconstruction constraint multipliers of any load channel, all updated candidate modes for that load channel are superimposed to obtain the current reconstructed phase sequence. The current reconstructed phase sequence is subtracted from the normalized phase sequence of that load channel to obtain the current signal reconstruction residual. The current signal reconstruction residual is multiplied by the multiplier update step size to obtain the multiplier correction. The multiplier correction is added to the reconstruction constraint multipliers of the previous round to obtain the current round of reconstruction constraint multipliers.

[0085] The multiplier update step size is set based on the change in signal reconstruction error during continuous iteration. Specifically, the step size is set starting from 1 for the candidate multiplier update. After completing the current round of multiplier update using the candidate multiplier update step size, the signal reconstruction error is recalculated. If the current round of signal reconstruction error is not greater than the previous round's error, the current candidate multiplier update step size is retained. If the current round of signal reconstruction error is greater than the previous round's error, the candidate multiplier update step size is divided by 2, and the current round of multiplier update is performed again. This process is repeated until the current round's signal reconstruction error is not greater than the previous round's error. For example, in this embodiment, the multiplier update step size is preferably between 0.1 and 1.

[0086] The number of candidate modes is set based on the number of rotating phase sampling points in the high-load section and the minimum number of phase sampling points required for a single mode to maintain stable oscillation. The minimum number of phase sampling points required for a single mode to maintain stable oscillation is set based on the shared center order and the estimation error of the modal phase. The specific setting method is as follows: On historical complete load data, starting with an oscillation cycle including 4 phase sampling points, the number of phase sampling points within one oscillation cycle is increased successively, and multivariate variational mode decomposition is performed for each cycle. The shared center order corresponding to the current number of phase sampling points is subtracted from the shared center order corresponding to the high-resolution reference result. The absolute value of the difference is taken and then divided by the reference shared center order to obtain the relative error of the shared center order. The modal phase relative error is calculated using the same method based on the current modal phase and the reference modal phase. When both the relative error of the shared center order and the relative error of the modal phase do not exceed the calibration relative error of the corresponding load channel, the number of phase sampling points that first meet the conditions is determined as the minimum number of phase sampling points.

[0087] The number of rotating phase sampling points in the high-load segment is divided by the minimum number of phase sampling points, and the quotient is rounded down to obtain the upper limit of the number of candidate modes. A number of candidate modes is generated between 2 and the upper limit. When the upper limit of the number of candidate modes exceeds the maximum number of modes that can be stably separated without generating repetitive modes in historical load data, the maximum number of modes is determined as the upper limit of the number of candidate modes. For example, in this embodiment, the number of candidate modes is preferably 2 to 8.

[0088] The bandwidth penalty parameter is set based on the modal reconstruction error, the degree of separation between modes, and the order stability of the shared center. Specifically, candidate bandwidth penalty parameters are generated at logarithmic intervals between 100 and 10000. Multivariate variational mode decomposition is then performed based on the number of candidate modes and each candidate bandwidth penalty parameter.

[0089] For any combination of parameters, sum all modes of the same load channel to obtain the reconstructed phase sequence. Subtract the original normalized phase sequence from the reconstructed phase sequence, square the differences, and sum them. Divide the sum of squared differences by the sum of squared values ​​of the original normalized phase sequence to obtain the reconstruction error.

[0090] Select any two modes, multiply and sum the values ​​of the two modes at the same rotational phase position to obtain the modal inner product. Square the values ​​of the two modes and sum them, then take the square root of each square to obtain the two modal norms. Divide the absolute value of the modal inner product by the product of the two modal norms to obtain the normalized inner product of the two modes. Summate the normalized inner products of all mode combinations and divide by the number of mode combinations to obtain the average normalized inner product.

[0091] The shared center order is updated based on the modal energy of the first and second halves of the high-load segment, respectively. The larger of the two shared center orders is subtracted from the smaller, and then divided by the average of the two shared center orders to obtain the percentage change in the shared center order. The maximum value among the reconstruction error, the average of the normalized inner product, and the percentage change in the shared center order is taken as the decomposition evaluation value. The candidate mode number and candidate bandwidth penalty parameter with the smallest decomposition evaluation value are selected as the target mode number and target bandwidth penalty parameters for the current high-load segment.

[0092] After each iteration of the alternating direction multiplier method, the relative change rate of all candidate modes is determined. The specific calculation process is as follows: Subtract the corresponding values ​​from the previous round from the values ​​of each candidate mode at each load channel and rotation phase position in the current round. Square the resulting differences and sum them to obtain the sum of squares of the modal update differences. Square and sum the values ​​of all candidate modes from the previous round to obtain the sum of squares of the modes from the previous round. Divide the sum of squares of the modal update differences by the sum of the sum of squares of the modes from the previous round and the lower limit of positive numbers to obtain the relative rate of change of the modes.

[0093] The positive lower limit is set based on the minimum non-zero square value of the standardized phase sequence, specifically: the minimum non-zero square value is determined as the positive lower limit. The convergence threshold is set based on the calibration relative error of each load channel, specifically taking the minimum value among the squared values ​​of the calibration relative errors of each load channel. Iteration stops when the modal relative rate of change is not greater than the convergence threshold for two consecutive rounds.

[0094] The upper limit of the number of iterations is set based on the maximum number of convergence rounds in the validation data. The specific method is as follows: Count the number of convergence rounds for all candidate parameter combinations in the validation data and extract the maximum number of convergence rounds. Divide the maximum number of convergence rounds by 10 and round the quotient up to obtain the reserved number of rounds. If the reserved number of rounds is less than 10, set 10 rounds as the reserved number of rounds. Add the maximum number of convergence rounds and the reserved number of rounds to obtain the upper limit of the number of iterations. Stop iterating when the alternating direction multiplier method reaches the upper limit of the number of iterations.

[0095] After iteration, inverse discrete Fourier transform is performed on each candidate mode frequency domain sequence to obtain the tool total torque mode phase sequence, axial load mode phase sequence, and drive current mode phase sequence.

[0096] The modal energy of the high-load segment is obtained by squaring and summing the modal values ​​within the high-load segment. The modal energy of the complete analysis interval is then obtained by squaring and summing the modal values ​​within the complete analysis interval. The high-load segment modal energy is then divided by the complete analysis interval modal energy to obtain the high-load energy percentage.

[0097] The high-load mode retention threshold is set based on the allowable false retention ratio of normal rock-breaking modes. Specifically, the setting method is as follows: divide the allowable number of falsely retained modes by the total number of normal rock-breaking modes to obtain the false retention ratio. Arrange the high-load energy proportions of normal rock-breaking modes within the in-phase reference window in ascending order. Extract the high-load energy proportion corresponding to when the cumulative proportion reaches 1 minus the false retention ratio to obtain the high-load mode retention threshold.

[0098] Modes with a high load energy percentage not lower than the high load mode retention threshold are retained. Based on the process stage division results, the rotational phase positions corresponding to the retained modes are associated with the load rise process, peak action process, and post-peak response process, respectively. Mode data belonging to the same load channel, the same mode number, and the same process stage, and with continuous rotational phase positions, are organized into the same high load mode segment to obtain the high load mode segment.

[0099] Step S303: Based on the high-load modal fragments, perform cross-channel homologous aggregation to obtain a set of high-load homologous fragments.

[0100] The specific steps of step S303 are as follows: Step S3031: Based on the high-load modal segments, extract the central order and modal phase, determine the signal channel identifier corresponding to each high-load modal segment, and associate the rotation phase range to obtain the modal segment feature set.

[0101] In this embodiment, the shared central order obtained from the multivariate variational mode decomposition is read as the central order of each high-load mode segment.

[0102] A Hilbert transform is performed on the high-load mode segment to obtain the analytic signal. The mode phase is determined based on the real and imaginary parts of the analytic signal. The specific calculation process is as follows: using the imaginary part as the ordinate and the real part as the abscissa, an arctangent operation is performed in four quadrants to obtain the initial mode phase at the current rotation phase position.

[0103] If the phase difference between the current initial modal phase and the previous rotational phase position is greater than half a cycle phase angle, subtract one full cycle phase angle from the current initial modal phase. When the phase difference is less than the negative half cycle phase angle, increase the current initial modal phase by one full cycle phase angle to obtain a continuous modal phase.

[0104] The source channel of the high-load modal segment is used to determine the tool total torque channel identifier, axial load channel identifier, or drive current channel identifier. The high-load segment identifier, process stage, center order, modal phase sequence, signal channel identifier, rotational phase range, and modal energy are correlated to obtain the modal segment feature set.

[0105] Step S3032: Based on the modal segment feature set, high-load modal segments with the same center order, overlapping rotation phase ranges, and different signal channel identifiers are combined to obtain cross-channel candidate segment groups.

[0106] In this embodiment, high-load modal segments are aggregated according to high-load segment identifiers and process stages.

[0107] The center-order consistency tolerance is set based on the center-order estimation error of the same shared mode in different load channels. The specific setting method is as follows: Based on the load data within the in-phase reference window, the center-order differences of the same shared mode in different load channels are statistically analyzed. All center-order differences are summed and then divided by the number of differences to obtain the average center-order difference. Each center-order difference is subtracted from the average, and the resulting difference is squared. All squared results are summed and divided by the number of differences minus 1. The square root of the result is then taken to obtain the standard deviation of the center-order difference. The average center-order difference is added to three times the standard deviation of the center-order difference to obtain the center-order consistency tolerance.

[0108] The difference in center orders is obtained by subtracting the smaller center order from the larger center order of the two modal segments. When the difference in center orders is not greater than the center order consistency tolerance, the two modal segments are considered to have the same center order.

[0109] Extract the larger start position and the smaller end position from the two rotation phase ranges. When the larger start position is not greater than the smaller end position, the range from the larger start position to the smaller end position is defined as the overlapping rotation phase range.

[0110] The minimum number of overlapping oscillation periods is set based on the convergence of the response phase difference fluctuation. Specifically, in historical stable segments of the same origin, the overlap range is increased periodically, starting from one complete oscillation period, and the response phase difference fluctuation is calculated for each period. The fluctuation after increasing the oscillation period is subtracted from the fluctuation before the increase, and the absolute value of the difference is taken. When the difference does not exceed the minimum modal phase resolution, the number of oscillation periods that first meet the condition is determined as the minimum number of overlapping oscillation periods.

[0111] The rotational phase length corresponding to a complete oscillation cycle is determined based on the average center order of the two modal segments. The rotational phase length is divided by the rotational phase division step size, and the result is rounded up. Then, it is multiplied by the minimum number of overlapping oscillation cycles to obtain the minimum number of overlapping sampling points.

[0112] High-load mode segments belonging to the same high-load section and the same process stage, with the same center order, overlapping rotation phase range, sufficient number of overlapping sampling points, and different signal channel identifiers are combined to obtain cross-channel candidate segment groups.

[0113] Step S3033: Based on the cross-channel candidate fragment group, calculate the response phase difference of each cross-channel candidate fragment within the overlapping rotation phase range, and determine the degree of fluctuation of the response phase difference.

[0114] In this embodiment, the modal phase sequences of different signal channels within the overlapping rotation phase range are aligned.

[0115] Rotate the phase position one by one, and subtract the modal phase of the second signal channel from the modal phase of the first signal channel to obtain the initial response phase difference. Subtract one full cycle phase angle from the initial response phase difference that is greater than half a cycle phase angle, and add one full cycle phase angle to the initial response phase difference that is less than the negative half cycle phase angle to obtain the response phase difference sequence.

[0116] The specific calculation process for the average response phase difference is as follows: Calculate the sine and cosine values ​​of each response phase difference separately. Add all the sine values ​​and divide by the number of response phase differences to obtain the sine average. Add all the cosine values ​​and divide by the number of response phase differences to obtain the cosine average. Using the sine average as the ordinate and the cosine average as the abscissa, perform an arctangent operation in four quadrants to obtain the average response phase difference.

[0117] The specific calculation process for the degree of response phase difference fluctuation is as follows: Calculate the minimum circumferential angular distance between each response phase difference and the average response phase difference. Square each minimum circumferential angular distance, add all the squared results, divide by the number of response phase differences, and take the square root of the result to obtain the degree of response phase difference fluctuation.

[0118] Step S3034: Based on the fluctuation of the response phase difference, determine the cross-stage stable relationship, and perform homogeneous fragment aggregation and attribute merging to obtain a high-load homogeneous fragment set.

[0119] In this embodiment, the cross-stage phase stability threshold is set based on the robust distribution of response phase difference fluctuations within the same central order range. Specifically, the fluctuations in response phase differences within the same central order range are arranged in ascending order, and the value corresponding to the middle position is extracted to obtain the median fluctuation. The median fluctuation is subtracted from each fluctuation, and the absolute value of the difference is taken. All absolute differences are then arranged in ascending order, and the absolute difference corresponding to the middle position is extracted to obtain the median absolute deviation. The median fluctuation and three times the median absolute deviation are added together to obtain the cross-stage phase stability threshold.

[0120] The minimum circumferential angular distance of the average response phase difference between adjacent process stages is calculated to obtain the cross-stage phase difference change. When the fluctuation of the response phase difference between two process stages and the cross-stage phase difference change are both no greater than the cross-stage phase stability threshold, a cross-stage stable relationship is determined to exist.

[0121] By using high-load modal segments as segment nodes, channel-based connections are established between segment nodes that meet the cross-channel stability condition, and stage-stable connections are established between segment nodes that meet the cross-stage stability condition, thus obtaining a channel-based segment connection structure.

[0122] Perform a connectivity traversal on the homologous segment connection structure, and divide segment nodes that can be reached from each other through channel homologous connections or stage-stable connections into the same connected component. Aggregate high-load modal segments within the same connected component into high-load homologous segments.

[0123] The specific calculation process for the merge center order is as follows: Multiply the center order of each modal segment by the corresponding modal energy to obtain the center order weight value. Sum all the center order weight values ​​and divide by the sum of all modal energies to obtain the merged center order. Take the union of the rotating phase range and the signal channel identifier to obtain the merged rotating phase range and the set of homologous channels, and then obtain the set of high-load homologous segments.

[0124] Step S304: Construct a phase correlation diagram of the high-load process based on the high-load homogeneous fragment set.

[0125] The specific steps of step S304 are as follows: Step S3041: Based on the set of high-load homologous fragments, map each high-load homologous fragment to a high-load action node, and associate it with the corresponding high-load section, rock breaking cycle, rotation phase range and process stage to obtain the set of high-load action nodes.

[0126] In this embodiment, a high-load homologous segment is mapped to a high-load action node. Node attributes include node identifier, high-load segment identifier, rock-breaking cycle, rotation phase range, process stage, central order, homologous channel set, and modal energy of each channel.

[0127] The representative value of the joint deviation is determined based on the median of the joint deviations within the corresponding rotation phase range. Specifically, the joint deviations are arranged in ascending order, and the joint deviation corresponding to the middle position is extracted. The specific calculation process for the change in propulsion displacement is as follows: the change in propulsion displacement is obtained by subtracting the propulsion displacement at the starting position from the propulsion displacement at the ending position of the rotation phase range.

[0128] The representative value of the joint deviation, the change in feed displacement, the tool rotation speed, and the feed speed are written into the node attributes to obtain the set of nodes under high load.

[0129] Step S3042: Based on the set of high-load action nodes, determine the node stage succession relationship within the same high-load segment, establish process evolution edges, and determine the high-load forming nodes, peak action nodes, and post-peak response nodes.

[0130] In this embodiment, nodes are grouped according to the high-load segment identifier, and arranged in the order of the load rise process, peak effect process, and post-peak response process. Within the same process stage, nodes are arranged according to the starting position of the rotation phase range.

[0131] Directed process evolution edges are established between adjacent nodes. The specific calculation process for the joint deviation change is as follows: subtract the representative value of the joint deviation of the preceding node from the representative value of the subsequent node's joint deviation to obtain the joint deviation change. The specific calculation process for the modal energy change is as follows: sum the modal energies of each channel of the subsequent and preceding nodes respectively, and then subtract the sum of the modal energies of the preceding node from the sum of the modal energies of the subsequent node to obtain the modal energy change. The specific calculation process for the node propulsion displacement difference is as follows: subtract the propulsion displacement change of the preceding node from the propulsion displacement change of the subsequent node to obtain the node propulsion displacement difference. The specific calculation process for the rotation phase interval is as follows: subtract the ending position of the rotation phase range of the preceding node from the starting position of the rotation phase range of the subsequent node to obtain the rotation phase interval. The joint deviation change, modal energy change, node propulsion displacement difference, and rotation phase interval are written into the process evolution edge attributes.

[0132] The node that first reaches the high load entry threshold in step S203 during the load increase process is identified as the high load forming node. The node with the largest representative value of the joint deviation during the peak action process is identified as the peak action node. The node containing the termination rotation phase of the post-peak response process is identified as the post-peak response node.

[0133] Step S3043: Based on the rock breaking cycle and rotation phase range of the post-peak response node, match the corresponding data from the in-phase data of adjacent rock breaking cycles, map them as in-phase verification nodes, and establish result verification edges.

[0134] In this embodiment, based on the rock breaking cycle and rotation phase range of the post-peak response node, load data, motion position data, tool rotation speed, feed speed and control execution data are extracted from the post-cycle in-phase data obtained in step S204 to obtain in-phase verification data.

[0135] The verification deviation sequence of the in-phase verification data is calculated using the method in step S202. The verification deviation sequence is then arranged in ascending order, and the verification deviation corresponding to the middle position is extracted to obtain the representative value of the verification deviation.

[0136] The change in propulsion displacement is obtained by subtracting the propulsion displacement at the starting position from the propulsion displacement at the termination position of the in-phase verification data. The in-phase verification data is then mapped to in-phase verification nodes.

[0137] Establish a result verification edge by using the post-peak response node as the starting node and the in-phase verification node as the ending node. Subtract the verification deviation representative value from the joint deviation representative value of the post-peak response node to obtain the in-phase load change. Subtract the propulsion displacement change of the post-peak response node from the verification propulsion displacement change to obtain the in-phase propulsion change, and write both types of changes into the result verification edge attributes.

[0138] Step S3044: Based on the control execution data, determine the execution response interval and the corresponding high-load action node, and establish the control action edge.

[0139] In this embodiment, the change in rotational speed command is obtained by subtracting the target tool rotational speed from the target rotational speed at the previous control moment from the target rotational speed at the next control moment. Similarly, the change in feed command is obtained by subtracting the target feed speed from the target feed speed at the previous control moment from the target feed speed at the next control moment. When the absolute value of any change is not less than the minimum command resolution value of the corresponding actuator, the next control moment is determined as the control command change moment.

[0140] The response initiation threshold is set based on the execution feedback fluctuations during the period when the control command remains unchanged. When calculating the average actual tool speed change, all actual tool speed changes during the period when the control command remains unchanged are summed, and then divided by the number of actual tool speed changes to obtain the average actual tool speed change. The average actual tool speed change is subtracted from each actual tool speed change to obtain the deviation of each speed change; each deviation is squared, and the sum of all squared deviations is divided by the number of actual tool speed changes minus 1 to obtain the mean square of the speed change deviations; the square root of the mean square of the speed change deviations is obtained to obtain the standard deviation of the actual tool speed change. Similarly, all actual feed rate changes are summed and divided by the number of actual feed rate changes to obtain the average actual feed rate change; the standard deviation of the actual feed rate change is obtained based on the difference between the actual feed rate change and the average actual feed rate change. The absolute value of the average actual tool speed change is added to three times the standard deviation of the actual tool speed change to obtain the speed response start threshold; the absolute value of the average actual feed rate change is added to three times the standard deviation of the actual feed rate change to obtain the feed response start threshold.

[0141] Search backwards from the moment the control command changes for sampling moments in which the actual change in tool rotation speed or actual feed rate continuously exceeds the corresponding response start threshold, and determine the first sampling moment that meets the condition as the response start moment.

[0142] The response stability determination threshold is set based on the tracking deviation distribution during the period when the control command remains unchanged. The specific setting method is as follows: calculate the median and median absolute deviation of the rotation tracking deviation and the propulsion tracking deviation respectively, and add the absolute value of the median to three times the median absolute deviation to obtain the corresponding response stability determination threshold.

[0143] Starting from the beginning of the execution response, search for sampling intervals where the tracking deviation is continuously no greater than the response stability threshold. When the number of consecutive stable sampling points reaches the final minimum number of consecutive sampling points determined in step S203, the first stable sampling time is determined as the execution response termination time, thus obtaining the execution response interval.

[0144] Extract high-load action nodes that overlap with the execution response interval. Use the nearest high-load action node before the control command change as the starting node of the control action edge, and the earliest high-load action node within the execution response interval as the ending node. Write the speed command change, feed command change, actual tool speed change, actual feed rate change, response delay, and tracking deviation into the control action edge attributes.

[0145] Step S3045: Based on the high-load action node set, in-phase verification nodes, process evolution edges, result verification edges, and control action edges, organize the graph structure to obtain the high-load process phase correlation graph.

[0146] In this embodiment, a node table is established based on high-load action nodes and in-phase verification nodes to record node type, high-load section, rock breaking cycle, rotation phase range, process stage, load characteristics, modal characteristics, motion position characteristics, and control execution characteristics.

[0147] An edge table is established based on process evolution edges, result verification edges, and control edges, recording edge type, starting node identifier, ending node identifier, and edge attributes. Outgoing edge indexes and incoming edge indexes are created for each node.

[0148] The high-load formation nodes, peak effect nodes, post-peak response nodes, and in-phase verification nodes belonging to the same high-load segment are organized into event subgraphs. All event subgraphs are arranged according to their occurrence time to obtain the high-load process phase correlation diagram.

[0149] Step S305: Based on the phase correlation diagram of the high-load process, perform recursive estimation of the switching state space under the closed lag constraint of the rock breaking result to obtain the high-load state result set.

[0150] The specific steps of step S305 are as follows: Step S3051: Based on the phase correlation diagram of the high-load process, extract the result verification path connecting the high-load forming node and the same-phase verification node.

[0151] In this embodiment, high-load forming nodes and in-phase verification nodes are extracted according to the high-load segment identifier. Starting from the high-load forming node, the process evolution edge is traversed to the post-peak response node, and then the result verification edge is traversed to the in-phase verification node to obtain the candidate result verification path.

[0152] The specific calculation process for stage completeness is as follows: sequentially check whether the candidate result verification path contains high load forming nodes, peak effect nodes, post-peak response nodes, and in-phase verification nodes. For each type of target node included, the stage completeness count is increased by 1 to obtain the stage completeness.

[0153] The specific calculation process for path phase continuity is as follows: Determine the larger starting phase and the smaller ending phase within the rotation phase range of adjacent nodes. When the larger starting phase is not greater than the smaller ending phase, subtract the larger starting phase from the smaller ending phase to obtain the phase overlap length of adjacent nodes. Add the phase overlap lengths of all adjacent nodes together and divide by the sum of the rotation phase range lengths of all nodes within the path to obtain the path phase continuity.

[0154] Select the candidate path with the highest stage completeness. When stage completeness is the same, select the candidate path with the highest path phase continuity to obtain the result verification path.

[0155] Step S3052: Take the high-load forming node as the state node to be confirmed, and combine it with the high-load homogeneous fragment set and control execution data to construct a switching state space model.

[0156] In this embodiment, the high load formation node is marked as a state node to be confirmed. A state to be confirmed indicates that a high load has been detected, but the post-peak response result, which characterizes whether effective crack propagation has occurred in the rock mass, and the in-phase verification result of adjacent rock-breaking cycles have not yet been fully formed. Existing state-space switching models typically update state probabilities based on current observation characteristics and directly use the state with the higher probability as the current state. However, both effective fracturing and ineffective obstruction states may exhibit a simultaneous increase in total tool torque, axial load, and drive current during the high load formation stage. It is difficult to accurately distinguish between the two states based solely on the current observation characteristics of the high load formation stage. To address this issue, this application sets the effective fracturing state and the ineffective obstruction state as two discrete states in the state-space switching model. Before the result verification path is closed, the high load formation node is kept as a state to be confirmed. The model first performs forward recursive state estimation based on the high load formation node, peak action node, and post-peak response node. After the in-phase verification node is formed, reverse state smoothing is performed based on the closed verification characteristics.

[0157] The state-space model includes discrete state variables, continuous state vectors, control input vectors, node observation vectors, continuous state transition relationships, node observation relationships, discrete state switching relationships, process noise parameters, and observation noise parameters. Discrete state variables include effective fracturing states and ineffective obstructed states. The continuous state vector comprises seven dimensions, representing, in order, the standardized joint deviation state, the standardized tool total torque mode energy state, the standardized axial load mode energy state, the standardized driving current mode energy state, the standardized shared center order state, the standardized response phase difference fluctuation state, and the standardized advance displacement change state. The control input vector comprises six dimensions, representing, in order, the rotational speed command change, the advance speed command change, the actual tool rotational speed change, the actual advance speed change, the rotational tracking deviation, and the advance tracking deviation. The node observation vector comprises seven dimensions, representing, in order, the standardized joint deviation representative value of the current node, the standardized tool total torque mode energy, the standardized axial load mode energy, the standardized driving current mode energy, the standardized shared center order, the standardized response phase difference fluctuation degree, and the standardized advance displacement change.

[0158] Based on the similar features in the closed-loop verification path, robust standardization is performed on the continuous state vector, control input vector, and node observation vector. For any feature dimension, all historical feature values ​​of that dimension are arranged in ascending order, and the feature value corresponding to the median position is extracted to obtain the median of that feature dimension. The median of that feature dimension is subtracted from each historical feature value, and the absolute value of the difference is taken. Then, all absolute differences are arranged in ascending order, and the absolute difference corresponding to the median position is extracted to obtain the median absolute deviation of that feature dimension. The median of that feature dimension is subtracted from the current feature value, and then divided by the median absolute deviation of that feature dimension to obtain the standardized feature value. When the median absolute deviation of that feature dimension is less than the minimum resolution value of the corresponding data, the minimum resolution value of the corresponding data is used as the denominator for standardization.

[0159] For both effective fracturing and ineffective hindered fracturing states, a 7x7 state transition matrix, a 7x6 control response matrix, a 7x7 observation mapping matrix, a 7x7 process noise covariance matrix, and a 7x7 observation noise covariance matrix are set up. During continuous state transitions, the 7-dimensional continuous state vector of the previous node is first multiplied by the state transition matrix corresponding to the current discrete state to obtain the continuous state prediction result without control. Then, the 6-dimensional control input vector of the current node is multiplied by the control response matrix corresponding to the current discrete state to obtain the control response result. The continuous state prediction result and the control response result without control are added according to their corresponding dimensions, and the process noise corresponding to the current discrete state is superimposed to obtain the continuous state of the current node. During node observations, the continuous state of the current node is multiplied by the corresponding observation mapping matrix to obtain the predicted node observation vector, and the observation noise corresponding to the current discrete state is superimposed to obtain the actual node observation vector.

[0160] Discrete state transitions are influenced by both the discrete state of the previous node and the control input vector of the current node. For the case where the previous node is in an effective fracturing state, two state transition scores are established: one for transitioning to an effective fracturing state and the other for transitioning to an ineffective, hindered state. For the case where the previous node is in an ineffective, hindered state, two state transition scores are established: one for transitioning to an effective fracturing state and the other for maintaining an ineffective, hindered state. Each state transition score is obtained by adding the corresponding state transition intercept to six control influence terms. Each control influence term is obtained by multiplying one dimension of the control input vector by the corresponding state transition coefficient. An exponential transformation is performed on the two state transition scores corresponding to the same previous state. Then, each exponential transformation result is divided by the sum of the two exponential transformation results corresponding to that previous state to obtain the probability of maintaining an effective fracturing state, the probability of transitioning from an effective fracturing state to an ineffective, hindered state, the probability of transitioning from an ineffective, hindered state to an effective fracturing state, and the probability of maintaining an ineffective, hindered state.

[0161] The initial values ​​for the continuous state are set based on the high-load forming nodes in the result verification path, and the 7-dimensional standardized node observation vectors of the high-load forming nodes are determined as the initial mean of the continuous state. The initial covariance of the continuous state is determined based on the 7-dimensional standardized node observation features of all high-load forming nodes in the closed result verification path. For any feature dimension, the average value of the corresponding feature values ​​of all high-load forming nodes is first calculated. Then, the average value is subtracted from each feature value, and the difference is squared. All squared results are summed, and then divided by the number of feature values ​​minus 1 to obtain the initial variance of that feature dimension. The initial variances of the seven feature dimensions are written to the main diagonal positions of the initial covariance matrix, and the remaining positions are set to zero, resulting in a 7-row, 7-column initial covariance matrix for the continuous state. When any initial variance is less than the square of the minimum resolution value of the corresponding standardized feature, the square of the minimum resolution value is used as the initial variance of that dimension.

[0162] The discrete states and model parameters are initialized based on the closed-loop result verification paths. From each closed-loop result verification path, the post-peak load release ratio, post-peak modal energy change ratio, adjacent-cycle in-phase load change ratio, and adjacent-cycle in-phase propagation change are extracted to obtain a 4-dimensional closed-loop initialization feature vector. Robust standardization of the 4-dimensional closed-loop initialization features yields a standardized closed-loop initialization feature vector. A K-means clustering algorithm with two cluster categories is used to cluster the standardized closed-loop initialization feature vector. The first initial cluster center is determined based on the upper quartiles of the four standardized closed-loop initialization features, and the second initial cluster center is determined based on the lower quartiles. Specifically, all standardized feature values ​​of each dimension are arranged in ascending order, and the values ​​corresponding to cumulative ratios reaching 0.75 and 0.25 are extracted. The four upper quartiles are arranged sequentially to obtain the first initial cluster center, and the four lower quartiles are arranged sequentially to obtain the second initial cluster center.

[0163] In any round of clustering iteration, for any closed result verification path, the standardized feature values ​​of its four dimensions are subtracted from the corresponding dimensions of the two cluster centers. The four differences are squared and summed to obtain the squared distance with the two cluster centers. The cluster center with the smaller squared distance is determined as the current category of the result verification path. When the squared distances with the two cluster centers are the same, the result verification path with the larger sum of the post-peak load release ratio and the ratio of in-phase load change in adjacent periods is assigned to category 1, and the remaining result verification paths are assigned to category 2. After completing the category assignment of all result verification paths, according to the feature dimensions, all the standardized feature values ​​currently assigned to the same category are added together and then divided by the number of result verification paths in that category to obtain the updated cluster centers.

[0164] The clustering stopping threshold is set based on the minimum effective change of the four standardized closed initialization features. Specifically, the minimum measured resolution of each original feature is divided by the robust standardized denominator for that dimension to obtain the standardized resolution. The minimum of the four standardized resolutions is then taken as the clustering stopping threshold. After updating the cluster centers, the corresponding values ​​from the previous round are subtracted from the current round's cluster center values ​​for each of the four dimensions. The four differences are squared and summed. The square root of the sum of the squares of the four differences is then taken to obtain the cluster center movement distance. Clustering stops when the movement distances of two cluster centers are not greater than the clustering stopping threshold, or when the categories of all result verification paths have not changed in two adjacent rounds.

[0165] The upper limit of the number of clustering iterations is set based on the maximum number of clustering convergence rounds of historical closed samples. The maximum number of clustering convergence rounds is divided by 10 and rounded up to obtain the number of reserved clustering rounds. Then, the maximum number of clustering convergence rounds and the number of reserved clustering rounds are added together to obtain the upper limit of the number of clustering iterations.

[0166] After clustering, the average values ​​of the post-peak load release ratio, post-peak modal energy change ratio, in-phase load change ratio in adjacent cycles, and in-phase propulsion change in adjacent cycles are calculated for each of the two categories. The category with relatively large average values ​​of the post-peak load release ratio, post-peak modal energy change ratio, and in-phase load change ratio in adjacent cycles, and with a positive average in-phase propulsion change in adjacent cycles, is identified as an effective fracturing state. The other category is identified as an ineffective hindered state. The initial probability of an effective fracturing state is obtained by dividing the number of result verification paths corresponding to the effective fracturing state by the total number of result verification paths; similarly, the initial probability of an ineffective hindered state is obtained by dividing the number of result verification paths corresponding to the ineffective hindered state by the total number of result verification paths.

[0167] Initial state labels obtained using the K-means clustering algorithm with a cluster size of 2 are used to initialize the state transition matrix, control response matrix, and observation mapping matrix. For any discrete state, all neighboring node samples initially labeled as that discrete state are extracted. The 7-dimensional continuous state vector of the previous node and the 6-dimensional control input vector of the current node are concatenated sequentially to obtain a 13-dimensional joint input vector. The 7-dimensional continuous state vector of the current node is determined as the corresponding output vector. Each joint input vector is multiplied by its transpose and accumulated to obtain the joint input accumulation matrix. Each current node continuous state vector is multiplied by its transpose and accumulated to obtain the state input correlation accumulation matrix. A regularization parameter is added to the main diagonal of the joint input accumulation matrix to obtain the regularized joint input matrix. The regularized joint input matrix is ​​inverted, and the state input correlation accumulation matrix is ​​multiplied by the inverse of the regularized joint input matrix to obtain a 7-row, 13-column combined parameter matrix. The first 7 columns are determined as the state transition matrix, and the last 6 columns are determined as the control response matrix. Using the same weighted least squares process, the observation mapping matrix is ​​obtained based on the current node's continuous state vector and the node's observation vector.

[0168] The regularization parameter is set based on the condition number of the matrix to be inverted. It is incremented incrementally from 0, and the parameter is added to the main diagonal of the matrix. The maximum and minimum singular values ​​of the regularized matrix are calculated, and the condition number is obtained by dividing the maximum by the minimum. The regularization parameter is stopped increasing when the condition number does not exceed the upper limit. The upper limit of the condition number is set based on the impact of the matrix inversion error on the observation and prediction results. The upper limit of the candidate condition number is incremented incrementally on the validation data, and the observation and prediction results are calculated based on the corresponding parameter matrices. The actual observation results are subtracted from the observation and prediction results, and the absolute value of the difference in each dimension is taken. The upper limit of the candidate condition number when the observation and prediction errors in each dimension do not exceed the corresponding feature calibration error is selected as the final upper limit of the condition number.

[0169] Based on the initialized state transition matrix and control response matrix, the continuous state prediction results for each node are calculated. The predicted continuous state results are subtracted from the actual continuous state results to obtain the process residual vector. Each process residual vector is then multiplied by its transpose and summed, and finally divided by the corresponding sample size to obtain the process noise covariance matrix for the corresponding discrete state. Based on the initialized observation mapping matrix, the observation vectors for the predicted nodes are calculated. The predicted node observation vectors are subtracted from the actual node observation vectors to obtain the observation residual vectors. Each observation residual vector is then multiplied by its transpose and summed, and finally divided by the corresponding sample size to obtain the observation noise covariance matrix for the corresponding discrete state.

[0170] Based on the initialized model parameters, the expectation-maximization algorithm is used for iterative updating of the model parameters. Each round of expectation-maximization iteration includes expectation processing and maximization processing. In expectation processing, forward recursive state estimation is performed sequentially for each closed result verification path. For any node and any discrete state, the continuous state estimate of the previous node is multiplied by the corresponding state transition matrix, and the control input vector of the current node is multiplied by the corresponding control response matrix. The results of the state transition matrix multiplication and the control response matrix multiplication are then added to obtain the continuous state prediction mean. The continuous state covariance of the previous node is multiplied sequentially by the state transition matrix and its transpose, and then added to the process noise covariance matrix to obtain the continuous state prediction covariance. The continuous state prediction mean is multiplied by the observation mapping matrix to obtain the predicted node observation vector. The predicted node observation vector is subtracted from the actual node observation vector to obtain the observation residual vector. The continuous state prediction covariance is multiplied sequentially by the observation mapping matrix and its transpose, and then added to the observation noise covariance matrix to obtain the observation residual covariance matrix.

[0171] The observation likelihood of the corresponding discrete state is calculated based on the observation residual vector and the observation residual covariance matrix. First, the observation residual covariance matrix is ​​inverted. Then, matrix multiplication is performed sequentially on the transpose of the observation residual vector, the inverse of the observation residual covariance matrix, and the observation residual vector to obtain the weighted distance of the observation residuals. This weighted distance is multiplied by -0.5 and subjected to an exponential transformation to obtain the observation residual exponent. The determinant of the observation residual covariance matrix is ​​calculated and its square root is taken to obtain the covariance scaling term. Based on the normalized scaling and covariance scaling term corresponding to the 7-dimensional Gaussian distribution, the observation residual exponent is normalized to obtain the observation likelihood.

[0172] Based on the current control input vector, four types of state transition probabilities are calculated. The probabilities of the two states from the previous node are multiplied by their corresponding state transition probabilities. The two products leading to the same current state are summed to obtain the predicted state probability for that current state. Each predicted state probability is multiplied by its corresponding observation likelihood to obtain two unnormalized forward state probabilities. These two unnormalized forward state probabilities are then divided by their sum to obtain the two types of forward state probabilities for the current node. Kalman update weights are calculated based on the continuous state prediction covariance and the observation residual covariance. The Kalman update weights are then multiplied by the observation residual vector to obtain the continuous state correction, which is then added to the continuous state prediction mean to obtain the continuous state update result for the corresponding discrete state.

[0173] After completing the forward recursion along the result verification path for all nodes, backward probability recursion is performed from the in-phase verification node to the high-load forming node. The two types of backward probabilities corresponding to the in-phase verification node are initialized to 1. For any preceding node and any preceding discrete state, the probability of transitioning from the preceding discrete state to the next node's two types of discrete states, the corresponding observation likelihood, and the corresponding backward probability are multiplied, and then the two products are added together to obtain the backward probability of the preceding discrete state. The two types of backward probabilities for the same node are then normalized. The forward state probability of any node is multiplied by its corresponding backward probability and normalized to obtain the two types of state smoothing probabilities for that node. The forward probability of the preceding node, the corresponding state transition probability, the observation likelihood of the next node, and the backward probability of the next node are multiplied and normalized to obtain the four types of state transition smoothing probabilities.

[0174] In the maximization process, the state smoothing probability is used as the sample weight to update the state transition matrix, control response matrix, observation mapping matrix, process noise covariance matrix, and observation noise covariance matrix. When updating the state transition matrix and control response matrix, the 13-dimensional joint input vector of each sample and the 7-dimensional continuous state vector of the current node are multiplied by the corresponding state smoothing probability, and then the combined parameter matrix is ​​recalculated according to the aforementioned weighted least squares process. When updating the observation mapping matrix, the state smoothing probability is used as the weight, and the observation mapping matrix is ​​recalculated based on the continuous state vector and the actual node observation vector. When updating the noise covariance matrix, each residual self-multiplication matrix is ​​multiplied by the corresponding state smoothing probability, and then all weighted residual self-multiplication matrices are summed and divided by the sum of all state smoothing probabilities.

[0175] Discrete state transition parameters are updated based on the state transition smoothing probabilities. For any previous discrete state, the state transition smoothing probabilities of transitioning to the two current discrete states are used as weights to establish a weighted binary logistic regression. The input is a 6-dimensional control input vector, and the parameters to be determined are two sets of state transition intercepts and corresponding control influence coefficients. Newton's iteration method is used to update the parameters. The two-class state transition probabilities of each sample are calculated based on the current parameters. The probability residual is obtained by subtracting the current state transition probability from the state transition smoothing probability. Combined with the control input vector and sample weights, the parameter update gradient and parameter update curvature matrix are obtained. The parameter update curvature matrix is ​​inverted and multiplied with the parameter update gradient to obtain the parameter correction amount. The state transition parameters are then updated based on the parameter correction amount. Newton's iteration stops when the absolute value of each dimension of the parameter correction amount is not greater than the state transition parameter stopping threshold.

[0176] After completing one expectation processing and one maximization processing, the log-likelihood of all closed result verification paths in the current round is calculated. For any node, the predicted probabilities of the two types of states are multiplied by the corresponding observation likelihoods to obtain the joint likelihood of the two types of nodes. The joint likelihoods of the two types of nodes are added together and transformed by the natural logarithm to obtain the node's log-likelihood. Then, the log-likelihoods of all nodes are added together to obtain the current round's model log-likelihood. The current round's model log-likelihood is subtracted from the previous round's model log-likelihood. The absolute value of the difference is taken and then divided by the sum of the absolute value of the previous round's model log-likelihood and the positive lower bound to obtain the relative rate of change of log-likelihood. The model stopping threshold is set based on the change in observed predictions caused by parameter changes in two adjacent rounds. Iteration stops when the relative rate of change of log-likelihood does not exceed the model stopping threshold for two consecutive rounds; iteration stops when the upper limit of the number of model iterations is reached, determined by the maximum number of convergence rounds and the reserved number of rounds based on the verification data, resulting in a switched state-space model.

[0177] For example, the switching state-space model uses a 7-dimensional continuous state vector, a 6-dimensional control input vector, and a 7-dimensional node observation vector, and sets two discrete states: an effective fracturing state and an ineffective obstructed state. The closed result verification path of Group A includes 240 records. The upper limit of the number of conditions determined based on the verification data is 100, the model stopping threshold is 0.00001, and the upper limit of the number of model iterations is 100 rounds. After the 37th iteration, the log-likelihood relative change rate is lower than the model stopping threshold for two consecutive rounds, and the model iteration is stopped.

[0178] Step S3053: Based on the result verification path and the switching state space model, perform forward state recursive estimation on the state node to be confirmed to obtain the forward state estimation result.

[0179] The specific steps of step S3053 are as follows: Step S30531: Based on the process evolution edge, organize the high load forming nodes, peak effect nodes and post-peak response nodes in sequence to obtain the node evolution sequence.

[0180] In this embodiment, high-load formation nodes, peak effect nodes, and post-peak response nodes are extracted along the result verification path and arranged according to the process evolution direction. When multiple nodes exist in the same process stage, they are arranged according to the starting position of the rotation phase range. When the starting positions of the rotation phase range are the same, they are arranged according to the sampling time to obtain the node evolution sequence.

[0181] Step S30532: Based on the node evolution sequence and control action edge, associate the control execution data with the corresponding high-load action node to obtain the node observation input sequence.

[0182] In this embodiment, the representative value of the joint deviation, modal energy, central order, response phase difference fluctuation degree and propulsion displacement change of each node are extracted sequentially to obtain the node observation feature sequence.

[0183] Based on the terminal node identifier of the control action edge, the changes in rotation speed command, feed command, actual tool rotation speed, actual feed speed, response delay, rotation tracking deviation, and feed tracking deviation are associated with the corresponding nodes to obtain the node control feature sequence.

[0184] When a node corresponds to multiple control edges, the control weight is determined based on the overlap length between the execution response interval and the node's time range. The specific calculation process is as follows: Determine the later start time and the earlier end time within both the execution response interval and the node's time range. Subtract the later start time from the earlier end time to obtain the overlap length. Divide the overlap length of a single control edge by the sum of the overlap lengths of all control edges to obtain the corresponding control weight. Multiply each control edge attribute by its corresponding control weight, and then sum the products of similar attributes to obtain the node's control characteristics.

[0185] When a node does not have a corresponding control action edge, read the most recent valid control command and executor feedback before the node occurred, and set the change in control command to zero.

[0186] The node observation feature sequence and the node control feature sequence are combined according to the node evolution order to obtain the node observation input sequence.

[0187] Step S30533: Based on the switching state space model, node evolution sequence and node observation input sequence, perform forward state recursive estimation to obtain node state probability sequence and state transition probability sequence.

[0188] In this embodiment, based on the switching state-space model obtained in step S3052, the state transition matrix, control response matrix, observation mapping matrix, process noise covariance matrix, observation noise covariance matrix, state transition intercept, and control influence coefficient corresponding to the effective fracturing state and the ineffective obstructed state are read, respectively. Based on the 7-dimensional standardized node observation vector of the high-load forming nodes, the initial mean of the continuous state is determined; based on the 7-row, 7-column initial covariance matrix of the continuous state determined in step S3052, the initial covariance of the continuous state is determined; based on the number of high-load forming nodes in the effective fracturing state and the number of high-load forming nodes in the ineffective obstructed state in the closed result verification path, the initial probability of the discrete state is determined. Specifically, the initial probability of the effective fracturing state is obtained by dividing the number of high-load forming nodes in the effective fracturing state by the total number of high-load forming nodes, and the initial probability of the ineffective obstructed state is obtained by dividing the number of high-load forming nodes in the ineffective obstructed state by the total number of high-load forming nodes.

[0189] Each node in the node evolution sequence is taken as the current node, and the corresponding 6-dimensional control input vector and 7-dimensional node observation vector are extracted. Based on the current node's control input vector, the scores for maintaining the effective fracturing state, transitioning from an effective fracturing state to an ineffective and obstructed state, transitioning from an ineffective and obstructed state to an effective fracturing state, and maintaining the ineffective and obstructed state are calculated respectively. When calculating any state transition score, the corresponding state transition intercept is added to the six control influence terms. Each control influence term is obtained by multiplying one dimension of the current control input vector by the corresponding control influence coefficient. The two state transition scores corresponding to the same previous state are subjected to exponential transformation, and each exponential transformation result is divided by the sum of the two exponential transformation results corresponding to the previous state to obtain the four types of state transition probabilities.

[0190] Based on the two types of state probabilities and four types of state transition probabilities of the previous node, the state prediction probability of the current node is calculated. The effective fracturing state probability of the previous node is multiplied by the effective fracturing state maintenance probability to obtain the first effective fracturing state prediction component; the ineffective obstructed state probability of the previous node is multiplied by the probability of the ineffective obstructed state transitioning to an effective fracturing state to obtain the second effective fracturing state prediction component; the two effective fracturing state prediction components are added together to obtain the current node's effective fracturing state prediction probability. Similarly, the effective fracturing state probability of the previous node is multiplied by the probability of the effective fracturing state transitioning to an ineffective obstructed state, and the ineffective obstructed state probability of the previous node is multiplied by the probability of the ineffective obstructed state maintenance probability; the two products are then added together to obtain the current node's ineffective obstructed state prediction probability.

[0191] Continuous state prediction is performed separately for effective fracturing and ineffective hindered fracturing states. For any discrete state, the continuous state estimate of the previous node is multiplied by the state transition matrix corresponding to the discrete state to obtain the state transition result; the control input vector of the current node is multiplied by the control response matrix corresponding to the discrete state to obtain the control response result; the two results are added according to their corresponding dimensions to obtain the mean of the continuous state prediction under the discrete state. The continuous state covariance of the previous node is multiplied sequentially by the state transition matrix and its transpose corresponding to the discrete state, and then added to the corresponding process noise covariance matrix to obtain the continuous state prediction covariance under the discrete state.

[0192] Based on the continuous state prediction mean and continuous state prediction covariance under any discrete state, the corresponding observation residuals and observation likelihoods are calculated. The continuous state prediction mean is multiplied by the corresponding observation mapping matrix to obtain the prediction node observation vector; the prediction node observation vector is subtracted from the actual node observation vector at the current node to obtain the observation residual vector; the continuous state prediction covariance is then multiplied sequentially by the observation mapping matrix and its transpose, and finally added to the corresponding observation noise covariance matrix to obtain the observation residual covariance matrix. Invert the observation residual covariance matrix, then perform matrix multiplication on the transpose of the observation residual vector, the inverse of the observation residual covariance matrix, and the observation residual vector to obtain the weighted distance of the observation residuals. Multiply the weighted distance of the observation residuals by -0.5 and perform an exponential transformation to obtain the observation residual exponent. Calculate the determinant of the observation residual covariance matrix and take its square root to obtain the covariance scaling term. Normalize the observation residual exponent based on the normalized scaling and covariance scaling term corresponding to the 7-dimensional Gaussian distribution to obtain the observation likelihood of the corresponding discrete state.

[0193] Multiplying the predicted probability of the effective fracturing state by the observed likelihood of the effective fracturing state yields the unnormalized probability of the effective fracturing state; multiplying the predicted probability of the ineffective obstructed state by the observed likelihood of the ineffective obstructed state yields the unnormalized probability of the ineffective obstructed state. Adding the two unnormalized probabilities gives the normalized denominator of the state probability. Dividing each of the two unnormalized probabilities by the normalized denominator of the state probability yields the current node's effective fracturing state probability and ineffective obstructed state probability. The lower bound of the probability denominator is set based on the smallest effective normalized denominator in the closed result verification path. Specifically, the smallest non-zero normalized denominator is divided by 10 to obtain the lower bound of the probability denominator. If the current state probability normalized denominator is less than the lower bound of the probability denominator, the current observed likelihood is not used to update the discrete state probability. Instead, the predicted probability of the effective fracturing state and the predicted probability of the ineffective obstructed state are divided by their sum to obtain the two state probabilities of the current node.

[0194] The continuous states corresponding to the two types of discrete states are updated separately. For any discrete state, the continuous state prediction covariance under that discrete state is multiplied by the transpose of the corresponding observation mapping matrix. The resulting matrix is ​​then multiplied by the inverse of the observation residual covariance matrix to obtain the Kalman update weights. The Kalman update weights are then multiplied by the observation residual vector to obtain the continuous state correction. The continuous state correction is then added to the continuous state prediction mean to obtain the continuous state update result under that discrete state. Finally, the Kalman update weights are multiplied by the observation mapping matrix, the resulting matrix is ​​subtracted from the 7x7 identity matrix, and then multiplied by the continuous state prediction covariance to obtain the continuous state update covariance under that discrete state.

[0195] The continuous state update result of the effective fracturing state is multiplied by the probability of the effective fracturing state at the current node, and the continuous state update result of the ineffective obstructed state is multiplied by the probability of the ineffective obstructed state at the current node. The two weighted results are then summed according to their corresponding dimensions to obtain the continuous state estimate of the current node. When calculating the continuous state covariance of the current node, the continuous state estimate is subtracted from the continuous state update results of the two types of discrete states to obtain two types of mean deviation vectors. Each mean deviation vector is multiplied by its transpose to obtain the corresponding mean deviation matrix. The continuous state update covariance and mean deviation matrix of each discrete state are summed and multiplied by the corresponding state probability to obtain the weighted result of the two types of covariance. The weighted result of the two types of covariance is then summed to obtain the continuous state covariance of the current node.

[0196] Based on the previous node's state probability, the four types of state transition probabilities, and the current node's observation likelihood, the joint probabilities of state transitions between the four adjacent nodes are calculated. The unnormalized joint probability of maintaining an effective fracturing state is obtained by multiplying the previous node's effective fracturing state probability, the effective fracturing state maintenance probability, and the current node's effective fracturing state observation likelihood. Similarly, the unnormalized joint probabilities corresponding to the transitions from an effective fracturing state to an ineffective / obstructed state, from an ineffective / obstructed state to an effective fracturing state, and maintaining an ineffective / obstructed state are obtained. The four types of unnormalized joint probabilities are summed to obtain the normalized denominator of the state transition joint probability. Then, the four types of unnormalized joint probabilities are divided by this normalized denominator to obtain the four types of state transition joint probabilities.

[0197] The probabilities of the current node's effective fracturing state and ineffective obstructed state are written into the node state probability sequence, and the joint probabilities of the four types of state transitions corresponding to the current node are written into the state transition probability sequence. The above forward state recursion process is repeated node by node along the node evolution sequence until the forward state recursion estimation of the post-peak response node is completed, thus obtaining the node state probability sequence and the state transition probability sequence.

[0198] Step S30534: Based on the node state probability sequence and state transition probability sequence, maintain a high load to form the unconfirmed state of the node before the result verification path is closed, and obtain the forward state estimation result.

[0199] In this embodiment, it is determined whether the in-phase verification node has been formed, and whether the corresponding load data, motion position data and control execution data have passed the validity check.

[0200] When a co-phase verification node has not yet been formed or the corresponding data contains invalid data intervals, the verification path for the determination result is not yet closed. The node state probability sequence and state transition probability sequence are retained, but the final state is not output based on the larger state probability. Nodes with high load formation are continued to be marked as pending confirmation.

[0201] The pending confirmation state is not a third rock mass state in the state-space model, but rather a setting of output constraints for state results for which closure evidence has not yet been obtained. The model still updates the probabilities of effective fracturing state and ineffective obstructed state.

[0202] When a co-phase verification node is formed and the corresponding data passes the validity check, the result verification path is considered closed. The result verification path, node state probability sequence, state transition probability sequence, and unconfirmed state identifier are correlated to obtain the forward state estimation result.

[0203] like Figures 2-4 As shown, step S3054: Based on the forward state estimation results, perform closed evidence-driven reverse state smoothing and control response correlation to obtain a high-load state result set.

[0204] The specific steps of step S3054 are as follows: Step S30541: Based on the forward state estimation results and result verification path, determine the state smoothing interval and organize the probability sequence to obtain the state smoothing interval and the probability sequence to be smoothed.

[0205] In this embodiment, the high-load forming node is determined as the starting node of the state smoothing interval, the in-phase verification node is determined as the ending node of the state smoothing interval, and all nodes passed along the result verification path between the two are included in the state smoothing interval.

[0206] Extract the effective fracturing state probability, ineffective obstructed state probability, and adjacent node state transition probability for each node within the state smoothing interval. Organize them in order from high-load forming nodes to in-phase verification nodes, and establish a reverse index from in-phase verification nodes to high-load forming nodes to obtain the probability sequence to be smoothed.

[0207] Step S30542: Based on the state smoothing interval, calculate the post-peak process response change and the in-phase load change of adjacent periods to obtain the closed-loop verification feature sequence.

[0208] In this embodiment, existing high-load identification methods typically calculate the load release ratio using peak load and post-peak load. Specifically, the peak load is subtracted from the post-peak load, and then divided by the peak load to obtain the load release ratio. This method can reflect the degree of load reduction after the peak, but tool deceleration, reduced feed speed, short-term actuator retraction, and partial tool detachment can all cause a decrease in post-peak load. Using the load release ratio alone cannot distinguish between temporary unloading and effective crack propagation.

[0209] To address this issue, this application uses post-peak load release, post-peak homogeneous mode energy change, in-phase load change in adjacent periods, and in-phase propagation change in adjacent periods as closure verification features.

[0210] The post-peak load release ratio is calculated based on the representative values ​​of the joint deviation of the peak-acting node and the post-peak response node. The specific calculation process is as follows: Subtract the representative value of the joint deviation of the post-peak response node from the representative value of the joint deviation of the peak-acting node to obtain the post-peak load release amount. Divide the post-peak load release amount by the representative value of the joint deviation of the peak-acting node to obtain the post-peak load release ratio.

[0211] Based on the modal energies of each load channel at the peak-acting node and the post-peak response node, the post-peak mode energy change ratio is calculated. The specific calculation process is as follows: Add the modal energies of each load channel at the peak-acting node to obtain the total peak-mode energy. Add the modal energies of each load channel at the post-peak response node to obtain the total post-peak mode energy. Subtract the total post-peak mode energy from the total peak-mode energy to obtain the post-peak mode energy change. Divide the post-peak mode energy change by the total peak-mode energy to obtain the post-peak mode energy change ratio.

[0212] Based on the representative value of the joint deviation of the post-peak response node and the in-phase verification node, the proportion of in-phase load change between adjacent periods is calculated. The specific calculation process is as follows: Subtract the representative value of the verification deviation of the in-phase verification node from the representative value of the joint deviation of the post-peak response node to obtain the in-phase load change between adjacent periods. Divide the in-phase load change between adjacent periods by the representative value of the joint deviation of the post-peak response node to obtain the proportion of in-phase load change between adjacent periods.

[0213] Based on the thrust displacement changes of the post-peak response node and the in-phase verification node, the in-phase thrust change in adjacent cycles is calculated. Specifically, the calculation process is as follows: subtract the thrust displacement change of the post-peak response node from the thrust displacement change of the in-phase verification node to obtain the in-phase thrust change in adjacent cycles.

[0214] The post-peak load release ratio, post-peak modal energy change ratio, adjacent cycle in-phase load change ratio, and adjacent cycle in-phase propulsion change are arranged according to the result verification path to obtain a closed verification feature sequence.

[0215] Step S30543: Based on the closure verification feature sequence, determine the valid fracturing closure evidence and the invalid blocked closure evidence to obtain the closure evidence results.

[0216] In this embodiment, existing technologies typically set separate thresholds for post-peak load release, propulsion displacement, and subsequent load change, and determine the state based on whether each feature exceeds the threshold. This method treats each feature as an independent hard determination condition, making it difficult to adapt to changes in feature distribution caused by different rock hardness, tool wear, and propulsion conditions.

[0217] To address this issue, this application calculates a weighted distance based on the variance of each closure verification feature in the two states, and then converts the weighted distance into soft evidence weights.

[0218] Based on the closed validation path, calculate the median and median absolute deviation of each closed validation feature. Subtract the corresponding median from the current feature value and then divide by the corresponding median absolute deviation to obtain the standardized closed validation feature.

[0219] The lower limit of the standardized denominator is set based on the measurement resolution value of the corresponding feature. Specifically, the measurement resolution value is determined as the lower limit of the denominator. When the median absolute deviation is less than the corresponding lower limit of the denominator, the lower limit of the denominator is used as the standardized denominator.

[0220] Based on the effective fracturing state observation parameters in the switched state-space model, the weighted distance of the effective fracturing state is calculated. The specific calculation process is as follows: Subtract the corresponding feature center of the effective fracturing state from each standardized closure verification feature to obtain the feature difference. Square each feature difference, then divide it by the variance of the corresponding feature under the effective fracturing state. Sum all the results to obtain the weighted distance of the effective fracturing state.

[0221] The lower limit of characteristic variance is set based on the square of the corresponding characteristic measurement resolution. When the state characteristic variance is less than the lower limit of characteristic variance, the lower limit of characteristic variance is used as the denominator for weighted distance calculation to avoid the denominator being too small, which would cause a single feature to be over-amplified.

[0222] Following the same method, the weighted distance of invalid and blocked states is obtained based on the feature center and corresponding feature variance of the invalid and blocked states.

[0223] Multiply the two weighted distances by -1 / 2 and then perform an exponential transformation to obtain the initial evidence strength for each of the two states. Divide the initial evidence strength for the effective fracturing state by the sum of the two initial evidence strengths to obtain the evidence weight for effective fracturing closure. Divide the initial evidence strength for the ineffective obstructed closure state by the sum of the two initial evidence strengths to obtain the evidence weight for ineffective obstructed closure.

[0224] By correlating the weights of valid fracturing closure evidence, invalid blocked closure evidence, and closure verification feature sequences, the closure evidence results are obtained.

[0225] Step S30544: Based on the probability sequence to be smoothed and the closure evidence results, perform reverse state smoothing along the result verification path to obtain the closed state probability sequence and the closed state transition probability sequence.

[0226] In this embodiment, traditional fixed-lag smoothing uses the smoothed state probability of the next node to correct the previous node. However, the smoothed state probability of the next node still mainly comes from the original state space model, making it difficult to use the in-phase verification results of adjacent periods to substantially correct the early state.

[0227] To address this issue, this application first uses the closed-loop evidence results to correct the state probability of the in-phase verification node, and then uses the corrected state probability as the starting point for reverse smoothing to sequentially correct the post-peak response node, peak effect node, and high load formation node along the result verification path.

[0228] Multiply the forward effective fracturing state probability of the in-phase verification node by the effective fracturing closure evidence weight to obtain the effective fracturing state correction value. Multiply the forward ineffective blocked state probability by the ineffective blocked closure evidence weight to obtain the ineffective blocked state correction value. Add the two correction values ​​to obtain the corrected normalized denominator. Divide the two correction values ​​by the corrected normalized denominator to obtain the two types of closure state probabilities of the in-phase verification node.

[0229] Starting from the node preceding the in-phase verification node, traverse backwards along the result verification path. For the current node to be smoothed, read the forward state probability of the current node, the state transition probability from the current node to the next node, the forward prediction probability of the next node, and the closed state probability of the next node.

[0230] For the current node's effective fracturing state, the probability of the effective fracturing state transitioning to the next node's effective fracturing state is multiplied by the probability of the next node's effective fracturing closure state, and then divided by the next node's effective fracturing forward prediction probability to obtain the first reverse correction component. The probability of the effective fracturing state transitioning to the next node's ineffective blocked state is multiplied by the probability of the next node's ineffective blocked closure state, and then divided by the next node's ineffective blocked forward prediction probability to obtain the second reverse correction component.

[0231] When the forward prediction probability of a certain state of the next node is zero, the corresponding backward correction component is set to zero to avoid zero-value division.

[0232] The first and second reverse correction components are added together, and then multiplied by the effective fracturing forward state probability of the current node to obtain the unnormalized smooth probability of the effective fracturing state. The unnormalized smooth probability of the ineffective obstructed state of the current node is calculated in the same manner.

[0233] Add the two unnormalized smoothed probabilities to obtain the smoothed normalized denominator. Divide the two unnormalized smoothed probabilities by the smoothed normalized denominator to obtain the probabilities of the two types of closed states of the current node.

[0234] Based on the closed-state probabilities and corresponding state transition probabilities of the current node and the next node, the joint probability of four types of state transitions is calculated. The specific calculation process is as follows: multiply the closed-state probability of the current node, its corresponding state transition probability, and the closed-state probability of the next node to obtain one type of joint probability of state transition. The four types of joint probabilities of state transitions are obtained in the same way. Each joint probability of state transition is then divided by the sum of the four types of joint probabilities to obtain the closed-state transition probability.

[0235] The reverse state smoothing of all nodes is completed along the result verification path to obtain the closed state probability sequence and the closed state transition probability sequence.

[0236] Step S30545: Based on the closed-state probability sequence and the closed-state transition probability sequence, perform high-load state determination and control response correlation to obtain the high-load state result set.

[0237] In this embodiment, the probabilities of effective fracturing closure and ineffective blocked closure are extracted for each high-load formation node. The state with the higher probability is determined as the state of the corresponding high-load section, and the higher probability is determined as the state confidence level.

[0238] Based on the occurrence time of high-load sections, the high-load sections identified as having effective fracturing states are arranged to obtain an effective fracturing state sequence. The high-load sections identified as having ineffective or hindered states are arranged to obtain an ineffective or hindered state sequence.

[0239] Based on the closed-state transition probability sequence and control action edge, the previous state, next state, speed command change, feed command change, actual tool speed change and actual feed speed change of adjacent high-load sections are extracted to obtain the control state transition record.

[0240] The discretization step size for the speed command change and the propulsion command change is set based on the minimum command resolution value of the corresponding actuator and the historical control change distribution. The specific setting method is as follows: Arrange the historical control changes in ascending order, extract the upper and lower quartiles respectively, and subtract the lower quartile from the upper quartile to obtain the interquartile range. Divide the interquartile range by the cube root of the number of historical control records to obtain the statistical grouping step size. Take the larger value between the statistical grouping step size and the minimum command resolution value of the corresponding actuator to obtain the discretization step size.

[0241] The speed command change and propulsion command change are divided into intervals according to the discretization step size, and the control state transition records are assigned to the corresponding control combination intervals.

[0242] The specific calculation process for the state transition response probability under each control combination is as follows: Add the corresponding closed-loop state transition probabilities within the same control combination interval to obtain the sum of the corresponding state transition probabilities. Add the four types of state transition probabilities separately to obtain the sum of all closed-loop state transition probabilities. Divide the sum of all types of state transition probabilities by the sum of all closed-loop state transition probabilities to obtain the probability of maintaining an effective fracturing state, the probability of an effective fracturing state transitioning to an ineffective / obstructed state, the probability of an ineffective / obstructed state transitioning to an effective fracturing state, and the probability of maintaining an ineffective / obstructed state.

[0243] By correlating the control combination interval, the probability of the four types of state transition responses, and the corresponding execution feedback changes, the high-load state transition response relationship is obtained. Combining the effective fracturing state sequence, the ineffective blocked state sequence, and the high-load state transition response relationship yields the high-load state result set.

[0244] exist Figure 2 The diagram in the upper left corner, showing the cutting tool acting on the crack propagation region, represents the post-peak response process. The crack formed below the tool indicates crack propagation in the rock mass after high load, corresponding to complete closure evidence. The solid waveform connected to it represents the response signal extracted from the post-peak response process. The bar chart to the right of the solid waveform represents the feature vector generated based on the post-peak load release ratio and the post-peak modal energy change ratio. The height of the bars indicates the relative magnitude of each feature, with high confidence indicating that the feature vector has a high degree of evidence completeness. The diagram in the lower left corner, showing the cutting tool acting on the fractured region, represents the in-phase verification process of adjacent rock-breaking cycles. The fragments below the tool represent the subsequent rock mass response under the same rotation phase, corresponding to incomplete closure evidence that has not yet formed a complete state conclusion independently. The dashed waveform connected to it represents the in-phase verification response signal of adjacent cycles. The bar chart to the right of the dashed waveform represents the feature vector generated based on the in-phase load change ratio and the in-phase propagation change of adjacent cycles, with low confidence indicating that this feature vector needs to participate in the closure evidence determination together with the post-peak response feature vector. Each solid arrow represents the correspondence between the response signal extracted from the rock breaking response process and the corresponding features calculated and organized from the response signal to verify the closure.

[0245] exist Figure 3 In the diagram, the upper rounded box represents valid fracturing closure evidence formed based on complete closure verification features. The single-peak curve within the box indicates a high degree of consistency between the current closure verification feature and the center of the effective fracturing state feature. The lower rounded box represents invalid obstructed closure evidence formed based on incomplete or persistent obstructed responses. The multi-peak curve within the box indicates that high loads have not been fully released or that in-phase loads in adjacent cycles remain at a high level. The dashed lines to the right of the two rounded boxes point to the plus-signed circle in the center, indicating that valid fracturing closure evidence and invalid obstructed closure evidence are converted into corresponding closure evidence weights and then normalized and weighted fused. The bar chart to the right of the plus-signed circle represents the fused evidence result. Each bar represents the relative strength of each closure verification feature or state evidence after fusion. The overall conclusion indicates that the fused result is used to correct the forward state probability of in-phase verification nodes.

[0246] exist Figure 4In the top, from right to left, are arranged the in-phase verification node, post-peak response node, peak effect node, and high-load formation node. Four hollow circles represent different nodes in the result verification path, and the dashed arrows between the hollow circles from right to left indicate the reverse propagation of closure evidence along the result verification path. The bar charts below each node represent the forward state probability or node observation characteristics of the corresponding node. The solid arrows pointing upwards from the bars to the corresponding nodes indicate that the node observation results are associated with the corresponding state nodes. The dashed lines below the bars converge and point in the reverse smoothing direction, indicating that the peak response node, peak effect node, and high-load formation node are corrected from right to left, starting with the closure state probability of the in-phase verification node and combining the forward state probabilities of each node and the state transition probabilities of adjacent nodes. The rounded corner box in the middle represents the reverse smoothing result, and the broken line formed by connecting multiple dots inside the box represents the corrected closure state probability sequence of each node in the result verification path. The checkmark graphic on the lower left indicates that the high-load state is determined based on the probability of the closed state, resulting in an effective fracturing state or an ineffective blocked state; the adjustment graphic on the lower right indicates that the closed state transition probability is correlated with the control execution data to obtain the control response relationship and control output command.

[0247] Figure 2 It is used to extract response signals from the post-peak response process and the in-phase verification process of adjacent rock breaking cycles, and to calculate the closed verification feature sequence; Figure 3 based on Figure 2 The obtained closure verification feature sequences are used to form valid fracturing closure evidence and invalid blocked closure evidence, respectively. Evidence weight fusion is then performed to obtain the closure evidence result. Figure 4 Will Figure 3 The obtained closure evidence results are introduced into the in-phase verification node, and reverse state smoothing is performed along the result verification path to obtain the closure state probability sequence and the closure state transition probability sequence. Then, high-load state determination and control response association are performed. The three figures correspond to the closure verification feature acquisition, closure evidence determination and fusion, and closure evidence-driven reverse state smoothing and control response association processes, respectively, and together represent the complete processing relationship of step S3054.

[0248] Step S4: Based on the high load state result set, perform rolling prediction of the next rock breaking phase state transition and closed-loop correction of the control quantity to obtain the load adaptive control instruction set.

[0249] Step S401: Based on the effective fracturing state sequence, the ineffective obstructed state sequence, the tool rotation speed, the feed speed, and the control execution data, determine the current high load state and the current control state to obtain the current state information.

[0250] In this embodiment, the effective fracturing state sequence and the ineffective blocked state sequence are uniformly sorted according to the occurrence time of the high-load section. The state corresponding to the high-load section with the latest occurrence time and which has completed the result verification path closure is determined as the current high-load state. If the current high-load section has not yet completed the path closure, it is kept in the pending confirmation state and does not overwrite the most recently closed high-load state.

[0251] Based on tool rotation speed, feed rate, and control execution data, extract the current actual tool rotation speed, current actual feed rate, current target tool rotation speed, current target feed rate, rotation tracking deviation, feed tracking deviation, rotation execution status, and feed execution status.

[0252] The upper limits for rotational tracking stability and propulsion tracking stability are set based on the tracking deviation distribution during the period when the control command remains constant. The specific setting method is as follows: calculate the median tracking deviation and the median absolute deviation respectively, and add the absolute value of the median tracking deviation to three times the median absolute deviation to obtain the corresponding upper limit for tracking stability.

[0253] When both the rotation and propulsion execution states are normal, and neither type of tracking deviation exceeds the corresponding tracking stability upper limit, the execution state is determined to be stable. If either execution state is a fault state, the execution state is determined to be faulty. All other situations are determined to be tracking adjustment states.

[0254] The rotational execution response coefficient is determined based on the change in rotational speed command and the actual change in tool rotational speed under stable execution conditions. The specific calculation process is as follows: Divide the actual change in tool rotational speed by the corresponding change in rotational speed command to obtain the single rotational response ratio. Arrange all valid single rotational response ratios in ascending order, and extract the response ratio corresponding to the middle position to obtain the rotational execution response coefficient. Records with a rotational speed command change of zero are not included in the calculation.

[0255] The propulsion execution response coefficient is determined in the same way, based on the change in propulsion speed command and the actual change in propulsion speed.

[0256] The current high load state, current actual tool speed, current actual feed rate, current target tool speed, current target feed rate, current execution state, tracking deviation, and execution response coefficient are correlated to obtain the current state information.

[0257] Step S402: Based on the current state information, generate a set of candidate control combinations for tool rotation speed and feed rate.

[0258] In this embodiment, the tool rotation speed and feed rate are extracted when the device protection is not triggered and the tracking deviation does not exceed the corresponding tracking stability limit under stable execution conditions, thus obtaining historical stable operation data.

[0259] The stable speed range is set based on historical stable operating data and the rated operating range of the rotary actuator. Specifically, the method is as follows: extract the minimum and maximum tool speeds from the historical stable operating data to obtain the historical stable speed range; then, determine the stable speed range by finding the intersection of the historical stable speed range and the rated operating range of the rotary actuator. The stable feed speed range is determined in the same way.

[0260] The candidate step size for the tool rotation speed is set based on the minimum effective command change of the rotary actuator and the historical stable rotation speed change distribution. Specifically, the setting method is as follows: extract the rotation speed command change that causes the actual tool rotation speed change to continuously exceed the rotation speed response start determination threshold in step S3044. Take the absolute value of each rotation speed command change and extract the smallest non-zero value to obtain the minimum effective command change of the rotary actuator.

[0261] Arrange the historical stable speed changes in ascending order and extract the upper and lower quartiles. Subtract the lower quartile from the upper quartile to obtain the interquartile range of speed changes. Divide the interquartile range of speed changes by the cube root of the number of historical stable speed change records to obtain the speed statistical step size. Take the larger value between the minimum effective command change of the rotary actuator and the speed statistical step size to obtain the candidate step size of the tool speed.

[0262] The candidate step size for feed rate is set based on the minimum effective command change of the feed actuator and the historical stable feed rate change distribution. The specific setting method is the same as that for the candidate step size for tool speed.

[0263] The maximum single adjustment amount of the stable speed is set based on the upper quartile of the speed command change under stable execution conditions. The specific setting method is as follows: take the absolute value of the speed command change and arrange it in ascending order; extract the value corresponding to the upper quartile statistical position to obtain the maximum single adjustment amount of the stable speed. The maximum single adjustment amount of the stable propulsion speed is determined in the same way.

[0264] When the current execution state is a stable execution state, candidate adjustment values ​​are generated based on the corresponding candidate step size within the range from the negative stable maximum single adjustment value to the positive stable maximum single adjustment value. Specifically, starting from the negative stable maximum single adjustment value, the candidate step size is increased successively until the positive stable maximum single adjustment value is reached.

[0265] When the current execution state is tracking and adjusting, the stable maximum single adjustment amount is multiplied by the corresponding execution response coefficient to obtain the reduced maximum single adjustment amount, and then a candidate adjustment amount is generated based on the corresponding candidate step size. When the current execution state is faulty, the candidate adjustment amount is set to zero.

[0266] The current actual tool speed is added to each candidate speed adjustment value to obtain the candidate tool speed value. The current actual feed rate is added to each candidate feed rate adjustment value to obtain the candidate feed rate value. Candidate values ​​that are outside the stable operating range are deleted. The remaining candidate tool speed and candidate feed rate values ​​are then combined in pairs to obtain the candidate control combination set.

[0267] Step S403: Based on the current high load state, candidate control combination set and high load state transition response relationship, perform state transition rolling prediction for the next rock breaking phase to obtain candidate state prediction result set.

[0268] In this embodiment, the candidate change in rotational speed is obtained by subtracting the current actual rotational speed from the candidate rotational speed. Similarly, the candidate change in feed rate is obtained by subtracting the current actual feed rate from the candidate feed rate.

[0269] Based on the high-load state transition response relationship, the control combination interval corresponding to each candidate control combination is found. When a candidate control combination falls into an existing control combination interval, the corresponding four types of state transition probabilities are read.

[0270] When a candidate control combination does not fall within an existing control combination interval, local interpolation is performed within the 2D control variable space comprised of candidate changes in engine speed and propulsion speed. A single control combination interval center in the 2D control variable space can only represent one discrete control position, and two control combination interval centers can only represent one local change direction. At least three non-collinear control combination interval centers are needed to form a local triangular neighborhood to simultaneously characterize the state transition probability changes in both the direction of candidate changes in engine speed and the direction of candidate changes in propulsion speed. Therefore, the three control combination intervals closest to the candidate control combination are selected from the existing control combination intervals, and it is determined whether the centers of these three control combination intervals are collinear.

[0271] To determine whether the centers of the three control combination intervals are collinear, the center of the first control combination interval is connected to the centers of the second and third control combination intervals, respectively, resulting in two direction vectors. Based on the components of these two direction vectors in the directions of candidate changes in rotational speed and propulsion speed, the area of ​​the triangle formed by the centers of the three control combination intervals is calculated. When the area of ​​the triangle is greater than zero, the centers of the three control combination intervals are determined to be non-collinear, and interpolation is performed using these three control combination intervals. When the area of ​​the triangle is equal to zero, the farthest control combination interval is replaced with the next control combination interval in ascending order of distance, until the centers of the three non-collinear control combination intervals are obtained. Therefore, the number of interpolation intervals is determined to be three, rather than simply increasing by one based on the dimension of the control variable.

[0272] The specific calculation process for the control combination distance is as follows: Subtract the center value of the speed change in the control combination interval from the candidate speed change, and then divide by the candidate tool speed step size to obtain the standardized difference of the speed. Subtract the center value of the feed rate change in the control combination interval from the candidate feed rate change, and then divide by the candidate feed rate step size to obtain the standardized difference of the feed rate. Square both standardized differences, add the two squared results, and then take the square root of the sum to obtain the control combination distance.

[0273] The specific calculation process for the interpolation weights is as follows: Calculate the reciprocal of the distances for the three control combinations, add the three reciprocals together to obtain the sum of the reciprocals. Divide the reciprocal of each control combination distance by the sum of the reciprocals to obtain the corresponding interpolation weight. When the distance of a certain control combination is zero, the state transition probability of that control combination interval is directly used, and no further interpolation is performed.

[0274] The state transition probabilities of the same type in the three control combination intervals are multiplied by their corresponding interpolation weights, and then the three products are added together to obtain the state transition probabilities corresponding to the candidate control combination.

[0275] When the current high-load state is an effective fracturing state, the probability of maintaining the effective fracturing state is determined as the predicted probability of the effective fracturing state in the next rock breaking phase, and the probability of the effective fracturing state turning into an ineffective and obstructed state is determined as the predicted probability of the ineffective and obstructed state in the next rock breaking phase.

[0276] When the current high-load state is an ineffective obstructed state, the probability of converting the ineffective obstructed state into an effective fracturing state is determined as the predicted probability of the effective fracturing state in the next rock breaking phase, and the probability of maintaining the ineffective obstructed state is determined as the predicted probability of the ineffective obstructed state in the next rock breaking phase.

[0277] Multiply the candidate change in rotational speed by the rotational execution response coefficient to obtain the predicted actual change in tool rotational speed. Multiply the candidate change in feed rate by the feed execution response coefficient to obtain the predicted actual change in feed rate. Add the predicted actual change in tool rotational speed to the current actual tool rotational speed to obtain the predicted tool rotational speed. Add the predicted actual change in feed rate to the current actual feed rate to obtain the predicted feed rate.

[0278] By correlating the candidate control combination, the prediction probabilities of the two types of states, the predicted tool rotation speed, and the predicted feed speed, a set of candidate state prediction results is obtained.

[0279] Step S404: Based on the current control state and the candidate state prediction result set, determine the target control combination and perform closed-loop correction of the control quantity to obtain the load adaptive control instruction set.

[0280] In this embodiment, the state prediction tolerance is set based on the prediction error of the switching state-space model in the validation data. Specifically, the setting method is as follows: subtract the corresponding actual closed state label from the state prediction probability of each validation sample. If the actual closed state belongs to the predicted state, the actual closed state label is set to 1; otherwise, it is set to 0. The absolute value of the difference is taken to obtain the absolute error of the state prediction. All absolute errors of the state prediction are arranged in ascending order, and the absolute error corresponding to the middle position is extracted to obtain the state prediction tolerance.

[0281] When the current high-load state is considered an effective fracturing state, the predicted probabilities of the effective fracturing state for each candidate control combination are compared to obtain the maximum predicted probability of the effective fracturing state. The predicted probability difference is obtained by subtracting the predicted probability of the effective fracturing state for each candidate control combination from the maximum predicted probability. Candidate control combinations whose predicted probability difference is not greater than the state prediction tolerance are retained.

[0282] When the current high-load state is an ineffective obstructed state, the predicted probabilities of the ineffective obstructed state transitioning to an effective fracturing state for each candidate control combination are compared to obtain the maximum state transition probability. The state transition probability difference is obtained by subtracting the state transition probability of each candidate control combination from the maximum state transition probability. Candidate control combinations whose state transition probability difference is not greater than the state prediction tolerance are retained, and then the candidate control combination with the lowest probability of maintaining the ineffective obstructed state is selected from these.

[0283] The specific calculation process for controlling the degree of change is as follows: Subtract the current actual tool speed from the candidate tool speed, take the absolute value of the difference, and then divide it by the maximum single adjustment amount of the stable speed to obtain the speed change ratio. Subtract the current actual feed speed from the candidate feed speed, take the absolute value of the difference, and then divide it by the maximum single adjustment amount of the stable feed speed to obtain the feed speed change ratio. Add the speed change ratio and the feed speed change ratio to obtain the degree of control change.

[0284] When multiple candidate control combinations with the same or similar state prediction results exist, the candidate control combination with the smallest degree of control change is selected to obtain the target control combination. The candidate values ​​for tool speed and feed rate are read from the target control combination. The candidate value for tool speed in the target control combination is determined as the target tool speed, and the candidate value for feed rate in the target control combination is determined as the target feed rate.

[0285] Based on the combination identifier of the target control combination, the predicted tool speed and predicted feed speed associated with the target control combination are extracted from the candidate state prediction result set, respectively, to obtain the target combination predicted tool speed and the target combination predicted feed speed. The target tool speed is subtracted from the target combination predicted tool speed to obtain the speed prediction tracking difference. The speed prediction tracking difference is divided by the rotational execution response coefficient to obtain the speed closed-loop compensation amount. The target feed speed is subtracted from the target combination predicted feed speed to obtain the feed speed prediction tracking difference. The feed speed prediction tracking difference is divided by the feed execution response coefficient to obtain the feed speed closed-loop compensation amount.

[0286] The effective lower limit of the spin response is set based on the lower tail distribution of the spin execution response coefficients under historical stable execution conditions. Specifically, the method is as follows: divide the number of allowable erroneous compensation records by the total number of historical stable response records to obtain the lower tail ratio. Arrange the historical spin execution response coefficients in ascending order, and extract the spin execution response coefficients corresponding to when the cumulative ratio reaches the lower tail ratio to obtain the effective lower limit of the spin response. The effective lower limit of the propulsion response is determined in the same way.

[0287] If the current rotational execution response coefficient is lower than the effective lower limit of the rotational response, the rotational speed closed-loop compensation is set to zero. If the current propulsion execution response coefficient is lower than the effective lower limit of the propulsion response, the propulsion speed closed-loop compensation is set to zero.

[0288] The target tool rotation speed and the rotation speed closed-loop compensation amount are added together to obtain the corrected tool rotation speed command. The target feed speed and the feed speed closed-loop compensation amount are added together to obtain the corrected feed speed command. The two corrected commands are then limited to the stable operating range determined in step S402.

[0289] The corrected tool speed command is quantified. The specific calculation process is as follows: Subtract the current actual tool speed from the corrected tool speed command to obtain the speed correction difference. Divide the speed correction difference by the candidate tool speed step size, and round the quotient to obtain the number of speed steps. Multiply the number of speed steps by the candidate tool speed step size to obtain the quantized speed correction amount. Add the quantized speed correction amount to the current actual tool speed to obtain the final rotation control command.

[0290] In the same manner, the corrected propulsion speed command is quantified based on the candidate propulsion speed step size to obtain the final propulsion control command.

[0291] By associating the final rotation control command, the final propulsion control command, the target control combination, the next rock-breaking phase state prediction result, and the current high load state, a load adaptive control command set is obtained.

[0292] Example 2: See Figure 5This embodiment provides a biomimetic rock-breaking tool load adaptive control system, including: a data acquisition module, a phase organization module, a recursive estimation module, and a prediction correction module.

[0293] The data acquisition module is used to collect the operational dataset of the bionic rock-breaking tool during the continuous rock-breaking process.

[0294] The phase organization module performs phase organization of the high-load process based on the running dataset to obtain the high-load rock breaking process sequence.

[0295] The recursive estimation module constructs a phase correlation diagram of the high-load process based on the high-load rock breaking process sequence, and performs recursive estimation of the switching state space under the closed lag constraint of the rock breaking result to obtain the high-load state result set.

[0296] The prediction and correction module, based on the high load state result set, performs rolling prediction of the next rock breaking phase state transition and closed-loop correction of the control quantity to obtain the load adaptive control instruction set.

[0297] The specific functions of each module described above are explained in the relevant content of the method in Embodiment 1, and will not be repeated here.

[0298] In addition, the parts of the technical solutions provided in the embodiments of this application that are consistent with the implementation principles of the corresponding technical solutions in the prior art have not been described in detail, so as to avoid excessive elaboration.

[0299] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the invention. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A biomimetic rock-breaking tool load adaptive control method, characterized in that, include: The operation dataset of the biomimetic rock-breaking tool during the continuous rock-breaking process is collected. The operation dataset includes load data, motion position data, tool rotation speed, feed speed, and control execution data. Based on the running dataset, phase organization of the high-load process is performed to obtain the high-load rock breaking process sequence; Based on the high-load rock breaking process sequence, a phase correlation diagram of the high-load process is constructed, and a recursive estimation of the switching state space under the closed lag constraint of the rock breaking result is performed to obtain a high-load state result set, which includes an effective fracturing state sequence, an ineffective obstructed state sequence, and a high-load state transition response relationship. Based on the high-load state result set, the rolling prediction of the next rock-breaking phase state transition and the closed-loop correction of the control quantity are performed to obtain the load adaptive control instruction set.

2. The biomimetic rock-breaking tool load adaptive control method according to claim 1, characterized in that, The process of organizing the phases of the high-load process based on the running dataset to obtain the high-load rock-breaking process sequence includes: Based on motion position data, the continuous rock-breaking process is divided into rock-breaking cycles, and the rotation phase within each rock-breaking cycle is determined to obtain the rock-breaking phase sequence. Based on the rock-breaking phase sequence, the in-phase load benchmark of the load data under each rotation phase is determined, and the deviation sequence is calculated; Based on the deviation sequence, high-load sections are extracted, and the rock-breaking cycle and rotation phase range corresponding to each high-load section are determined. Based on the rock-breaking cycle and rotation phase range corresponding to each high-load section, in-phase data of adjacent rock-breaking cycles are extracted and correlated with each high-load section to obtain the high-load rock-breaking process sequence.

3. The biomimetic rock-breaking tool load adaptive control method according to claim 2, characterized in that, Based on the high-load rock-breaking process sequence, a phase correlation diagram of the high-load process is constructed, and a recursive estimation of the switching state space under the closed lag constraint of the rock-breaking results is performed to obtain a high-load state result set, including: Based on the high-load rock breaking process sequence, the deviation sequence change corresponding to each high-load section was determined, and the load rising process, peak action process and post-peak response process were divided to obtain the process stage division results. Based on the high-load rock breaking process sequence and process stage division results, phase synchronization mode extraction was performed, and process stage correlation was performed to obtain high-load mode fragments; Based on high-load modal fragments, cross-channel homology aggregation is performed to obtain a set of high-load homology fragments; Construct a phase correlation graph for high-load processes based on a set of high-load homogeneous fragments; Based on the phase correlation diagram of the high-load process, a recursive estimation of the switching state space under the closed lag constraint of the rock breaking result is performed to obtain the high-load state result set.

4. The biomimetic rock-breaking tool load adaptive control method according to claim 3, characterized in that, The method of performing cross-channel homology aggregation based on high-load modal fragments to obtain a set of high-load homology fragments includes: Based on the high-load modal segments, the central order and modal phase are extracted, the signal channel identifier corresponding to each high-load modal segment is determined, and the rotation phase range is associated to obtain the modal segment feature set; Based on the modal segment feature set, high-load modal segments with the same center order, overlapping rotation phase range, and different signal channel identifiers are combined to obtain cross-channel candidate segment groups; Based on the cross-channel candidate fragment group, the response phase difference of each cross-channel candidate fragment within the overlapping rotation phase range is calculated, and the degree of fluctuation of the response phase difference is determined. Based on the degree of fluctuation of the response phase difference, cross-stage stable relationships are determined, and homogeneous fragment aggregation and attribute merging are performed to obtain a high-load homogeneous fragment set.

5. The biomimetic rock-breaking tool load adaptive control method according to claim 3, characterized in that, The construction of a phase correlation graph for a high-load process based on a set of high-load homogeneous fragments includes: Based on the set of high-load homologous fragments, each high-load homologous fragment is mapped to a high-load action node, and associated with the corresponding high-load section, rock breaking cycle, rotation phase range and process stage to obtain the set of high-load action nodes. Based on the set of high-load action nodes, the node stage succession relationship within the same high-load segment is determined, process evolution edges are established, and high-load formation nodes, peak action nodes, and post-peak response nodes are identified. Based on the rock-breaking cycle and rotation phase range of the post-peak response node, corresponding data are matched from the in-phase data of adjacent rock-breaking cycles, mapped to in-phase verification nodes, and result verification edges are established. Based on the control execution data, determine the execution response range and the corresponding high-load action node, and establish the control action edge; Based on the high-load action node set, in-phase verification nodes, process evolution edges, result verification edges, and control action edges, a graph structure is organized to obtain a high-load process phase correlation graph.

6. The biomimetic rock-breaking tool load adaptive control method according to claim 5, characterized in that, The method involves recursively estimating the switching state space under the closed-loop hysteresis constraint of the rock breaking results based on the phase correlation diagram of the high-load process, resulting in a high-load state result set, including: Based on the phase correlation graph of the high-load process, the result verification path connecting the high-load forming node and the same-phase verification node is extracted; The high-load forming node is taken as the state node to be confirmed, and a switching state space model is constructed by combining the high-load homogeneous fragment set and control execution data. Based on the result verification path and the switching state space model, a forward state recursive estimation is performed on the state node to be confirmed to obtain the forward state estimation result. Based on the forward state estimation results, closed evidence-driven backward state smoothing and control response correlation are performed to obtain a high-load state result set.

7. The biomimetic rock-breaking tool load adaptive control method according to claim 6, characterized in that, The process, based on the result verification path and the switching state space model, involves forward recursive estimation of the state node to be confirmed, yielding forward state estimation results, including: Based on the process evolution edge, the high load forming nodes, peak effect nodes, and post-peak response nodes are sequentially organized to obtain the node evolution sequence; Based on the node evolution sequence and control action edge, the control execution data is associated with the corresponding high-load action node to obtain the node observation input sequence; Based on the switching state space model, node evolution sequence and node state input sequence, forward state recursive estimation is performed to obtain node state probability sequence and state transition probability sequence. Based on the node state probability sequence and the state transition probability sequence, the node is kept in a state of unconfirmation under high load before the result verification path is closed, thus obtaining the forward state estimation result.

8. The biomimetic rock-breaking tool load adaptive control method according to claim 6, characterized in that, Based on the forward state estimation results, closed-evidence-driven backward state smoothing and control response correlation are performed to obtain a high-load state result set, including: Based on the forward state estimation results and result verification path, the state smoothing interval is determined and the probability sequence is organized to obtain the state smoothing interval and the probability sequence to be smoothed. Based on the state smoothing interval, the post-peak process response change and the in-phase load change of adjacent periods are calculated to obtain the closed-loop verification feature sequence. Based on the closure verification feature sequence, valid fracturing closure evidence and invalid blocked closure evidence are identified, and the closure evidence results are obtained. Based on the probability sequence to be smoothed and the closure evidence results, reverse state smoothing is performed along the result verification path to obtain the closed state probability sequence and the closed state transition probability sequence. Based on the closed-state probability sequence and the closed-state transition probability sequence, high-load state determination and control response correlation are performed to obtain a high-load state result set.

9. The biomimetic rock-breaking tool load adaptive control method according to claim 1, characterized in that, Based on the high-load state result set, the next rock-breaking phase state transition rolling prediction and control quantity closed-loop correction are performed to obtain the load adaptive control instruction set, including: Based on the effective fracturing state sequence, the ineffective obstructed state sequence, the tool rotation speed, the feed speed, and the control execution data, the current high load state and the current control state are determined, and the current state information is obtained. Based on the current state information, a set of candidate control combinations for tool rotation speed and feed rate is generated; Based on the current high load state, candidate control combination set and high load state transition response relationship, the state transition rolling prediction of the next rock breaking phase is carried out to obtain the candidate state prediction result set; Based on the current control state and the candidate state prediction result set, the target control combination is determined, and the control quantity is closed-loop corrected to obtain the load adaptive control instruction set.

10. A biomimetic rock-breaking tool load adaptive control system, used to implement the biomimetic rock-breaking tool load adaptive control method according to any one of claims 1 to 9, characterized in that, include: The data acquisition module is used to collect the operational dataset of the bionic rock-breaking tool during the continuous rock-breaking process; The phase organization module performs phase organization of the high-load process based on the running dataset to obtain the high-load rock breaking process sequence; The recursive estimation module constructs a phase correlation diagram of the high-load process based on the high-load rock breaking process sequence, and performs recursive estimation of the switching state space under the closed lag constraint of the rock breaking result to obtain the high-load state result set. The prediction and correction module, based on the high load state result set, performs rolling prediction of the next rock breaking phase state transition and closed-loop correction of the control quantity to obtain the load adaptive control instruction set.

Images