Fault prediction and self-protection method and system for heave compensation of diving bell

Abstract

The invention provides a fault prediction and self-protection method and system for heave compensation of a diving bell, and the method comprises the following steps: obtaining the displacement of a hydraulic cylinder of a heave compensation device and the pressure data of an energy accumulator, constructing a hysteresis loop, converting the hysteresis loop into a feature field image, and analyzing and calculating a topological feature descriptor through geometric moment. By performing inverse frequency domain analysis on the pressure data, a friction texture signal is separated and an energy entropy value of the friction texture signal is calculated. Mapping the topological feature descriptor to a control variable space defined by a nonlinear cuspidal point mutation model, and calculating a geometric distance from a current state point to the boundary of a model bifurcation set; and solving an instability prediction index by combining the distance and the energy entropy. And when the index exceeds the limit, the system triggers variable damping control and gradually locks the hydraulic cylinder to realize self-protection. According to the method, the technical problem that the sudden instability failure of the system caused by multi-factor coupling is difficult to warn in advance by a traditional method is solved.

Description

Technical Field

[0001] This invention belongs to the field of fault prediction technology, specifically relating to a fault prediction and self-protection method and system for diving bell heave compensation. Background Technology

[0002] As a critical piece of equipment in deep-sea operations, the diving bell experiences severe heave during deployment and retrieval due to the influence of surface waves. To ensure the diving bell's depth stability at the underwater operating point and the safety of the cable system, an active or passive heave compensation device must be installed. The core of this device is a complex hydraulic-pneumatic coupling system consisting of hydraulic cylinders, accumulators, and control valve assemblies. Its goal is to counteract the heave of the mother ship in real time, ensuring the stability of the diving bell. Operating long-term in the high-pressure, high-dynamic, and highly corrosive marine environment, the device inevitably experiences wear, aging, and performance degradation in its internal seals, hydraulic valve cores, sensors, and other critical components.

[0003] To ensure reliable operation, existing monitoring systems typically monitor key parameters in hydraulic systems in real time, such as pressure, displacement, and oil temperature, and set corresponding alarm thresholds. When a parameter exceeds a preset safety range, the system issues an alarm. This threshold-based method is effective for detecting significant faults that have already occurred (such as sudden pressure drops or excessively high oil temperatures). However, many serious faults in heave compensation devices do not manifest as a slow, linear exceedance of parameters, but rather as a sudden system instability. For example, seals may suddenly rupture after long-term wear, or control valves may suddenly become stuck due to oil contamination during high-frequency operation. The commonality of these faults is that, just before the failure occurs, the system's macroscopic parameters may still fluctuate within the set normal threshold range, but the deep dynamic characteristics within the system have undergone a qualitative change. Traditional threshold monitoring methods cannot detect this process of the system approaching the instability critical point from a stable state. Summary of the Invention

[0004] This invention provides a fault prediction and self-protection method and system for diving bell heave compensation to solve the above-mentioned technical problems.

[0005] In a first aspect, the present invention provides a fault prediction and self-protection method for diving bell heave compensation, the method comprising the following steps: Acquire real-time operating data of the heave compensation device, extract hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle from the real-time operating data, map them to a two-dimensional coordinate system to construct a closed hysteresis loop, and convert the closed hysteresis loop into a binary hysteresis feature field image. Geometric moment analysis is performed on the hysteresis feature field image, and topological feature descriptors for quantizing the discreteness and skewness characteristics of the hysteresis loop are calculated. Homomorphic deconvolution is used to process the accumulator pressure data and transform it to the cepstrum domain. Low-frequency components of the accumulator pressure data are filtered out by inverted high-pass filtering. At the same time, the friction texture signal characterizing the mechanical contact state is separated and the energy entropy value of the friction texture signal is calculated. The discreteness characteristics are mapped to splitting factor control variables, and the skewness characteristics are mapped to normal factor control variables. A control variable space consisting of splitting factor control variables and normal factor control variables is constructed, and the current state point of the heave compensation device in the control variable space is determined. The geometric distance from the current state point to the boundary of the bifurcation set of the nonlinear cusp catastrophe model is calculated based on the coordinates of the current state point in the control variable space. The instability prediction index is then calculated by combining the geometric distance and the energy entropy value. When the instability prediction index meets the preset risk threshold condition, a variable damping control command is generated to adjust the hydraulic throttling element of the heave compensation device, and the heave compensation hydraulic cylinder of the heave compensation device is gradually locked to complete the fault self-protection.

[0006] Optionally, the step of extracting hydraulic cylinder displacement data and accumulator pressure data from real-time operating data within a single heave cycle, mapping them to a two-dimensional coordinate system to construct a closed hysteresis loop, and converting the closed hysteresis loop into a binary hysteresis feature field image includes the following steps: Set a two-dimensional grid resolution that is adapted to the displacement range of the hydraulic cylinder and the pressure variation range of the accumulator, and initialize the two-dimensional grid. The two-dimensional grid resolution is used to define the pixel density of the hysteresis feature field image. The hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle in the real-time operation data are processed by cubic spline interpolation to generate a smooth and continuous closed hysteresis loop trajectory. The coordinates of the closed hysteresis loop trajectory are normalized to a unit interval, and the normalized coordinates are mapped to a two-dimensional grid. The grid cells passed by the closed hysteresis loop trajectory are marked as foreground pixels, and the grid cells not passed by the closed hysteresis loop trajectory are marked as background pixels, thus generating an initial two-dimensional grid image. The initial two-dimensional mesh image is preprocessed using morphological closing operations to obtain a two-dimensional mesh image. The two-dimensional mesh image is then transformed into a binary matrix containing only 0 and 1 values ​​to generate a hysteresis feature field image for feature extraction.

[0007] Optionally, the step of performing geometric moment analysis on the hysteresis feature field image and calculating the topological feature descriptor for quantizing the discreteness and skewness features of the hysteresis loop includes the following steps: The zeroth moment of the hysteresis feature field image is calculated to obtain the closed area of ​​the hysteresis loop, and the geometric centroid coordinates of the hysteresis loop in the two-dimensional coordinate system are determined by calculating the first moment of the hysteresis feature field image. The second and third central moments of the hysteresis feature field image are calculated based on the geometric centroid coordinates. The second and third central moments are then normalized using the zero-order moment to obtain the normalized central moments. The first invariant characteristic formula is constructed using the second-order central moment in the normalized central moments, and the sum of the moments of inertia in the orthogonal directions is defined as the discreteness characteristic that characterizes the discreteness of the figure. The second invariant feature formula is constructed using the third central moment in the normalized central moments. The sum of the projection skewness in the orthogonal directions is defined as the skewness feature that characterizes the degree of skewness of the figure. The calculated discreteness feature and the skewness feature are combined to form a topological feature descriptor.

[0008] Optionally, the step of processing the accumulator pressure data using homomorphic deconvolution and transforming the accumulator pressure data to the cepstrum domain, filtering out the low-frequency components of the accumulator pressure data through inverted high-pass filtering, and simultaneously separating the friction texture signal characterizing the mechanical contact state and calculating the energy entropy value of the friction texture signal includes the following steps: A Hamming window is applied to the collected accumulator pressure data, and a fast Fourier transform is performed on the windowed accumulator pressure data to obtain the pressure spectrum. Taking the natural logarithm of the amplitude of the pressure spectrum transforms the convolutional relationship in the time domain of the pressure spectrum into an additive relationship in the logarithmic amplitude spectrum in the frequency domain. Perform an inverse fast Fourier transform on the logarithmic amplitude spectrum in the frequency domain to map the accumulator pressure data to the inverted frequency domain to obtain the cepstrum sequence of the accumulator pressure data; The inverted frequency cutoff threshold is set according to the wave compensation period of the heave compensation device, and an inverted high-pass filter is constructed to filter out the low inverted frequency components representing wave motion and fluid pressure in the cepstral sequence. An inverse homomorphic transform is performed on the high inverse frequency residual component after processing by the inverse high-pass filter, and the high inverse frequency residual component is restored to the time domain to obtain the friction texture signal; The probability density function of the friction texture signal is calculated, and the entropy value of the probability density function is calculated based on the Shannon entropy formula to obtain the energy entropy value that quantifies the complexity of the mechanical contact surface.

[0009] Optionally, the step of setting the inverted frequency cutoff threshold according to the wave compensation period of the heave compensation device and constructing an inverted high-pass filter to filter out the low inverted frequency components representing wave motion and fluid pressure in the cepstral sequence includes the following steps: Calculate the amplitude envelope of the cepstrum sequence and normalize the amplitude envelope to construct a histogram of the cepstrum energy probability distribution; The information complexity of the distribution histogram is quantified by Raney entropy, and the sum of the Raney entropy of the background class cepstral distribution and the foreground class cepstral distribution is defined as the objective function of the candidate cepstral cutoff threshold. Traverse all candidate cepstral frequency cutoff thresholds within the domain of the cepstral sequence, solve the global maximum value of the objective function using a genetic algorithm, and determine the candidate cepstral frequency cutoff threshold corresponding to the global maximum value as the optimal cepstral frequency cutoff threshold; Based on the optimal cepstral frequency cutoff threshold, a cosine tapered window function with smooth transition characteristics is constructed. The cosine tapered window function is multiplied with the cepstral sequence to suppress wave fluid interference at low cepstral frequencies while preserving the friction characteristics of high cepstral frequencies.

[0010] Optionally, the step of mapping the discreteness characteristics to splitting factor control variables, mapping the skewness characteristics to normality factor control variables, constructing a control variable space composed of splitting factor control variables and normality factor control variables, and determining the current state point of the heave compensation device in the control variable space includes the following steps: A nonlinear cusp mutation model of the heave compensation device is constructed and a standard cusp mutation potential function is established. The standard cusp mutation potential function includes a state variable, a quartic term, a quadratic term weighted by a splitting factor, and a linear term weighted by a normality factor. Historical operating data of the heave compensation device under healthy conditions were collected, and a baseline distribution model was established by calculating the dispersion and skewness characteristics of the historical operating data. Calculate the first standardized Z-score of the discreteness feature extracted at the current time relative to the baseline distribution model, and the second standardized Z-score of the skewness feature extracted at the current time relative to the baseline distribution model; The first standardized Z-score is compressed and mapped using a nonlinear Sigmoid function to generate a splitting factor control variable corresponding to the quadratic term of the standard cusp mutation potential function. The hyperbolic tangent function is used to compress and map the second standardized Z-score, generating a normal factor control variable corresponding to the first term of the standard cusp mutation potential energy function; A control variable space consisting of splitting factor control variables and normal factor control variables is constructed, and the current state point of the heave compensation device in the control variable space is determined based on the standard cusp mutation potential energy function.

[0011] Optionally, the step of calculating the geometric distance from the current state point to the boundary of the bifurcation set of the nonlinear cusp catastrophe model based on the coordinates of the current state point in the control variable space, and then combining the geometric distance and the energy entropy value to calculate the instability prediction index includes the following steps: Based on the mathematical definition of the nonlinear cusp catastrophe model, a bifurcation set discrimination equation consisting of splitting factor control variables and normal factor control variables is constructed, and bifurcation set curves are plotted in the control variable space. An iterative search algorithm is used to calculate the shortest Euclidean distance from the coordinates of the current state point in the control variable space to the bifurcation set curve. After normalizing the energy entropy value, a weighting coefficient is assigned to the energy entropy value. The instability prediction index is calculated by combining the shortest Euclidean distance with the weighted energy entropy value and using the instability prediction index calculation formula.

[0012] Optionally, when the instability prediction index meets the preset risk threshold condition, generating a variable damping control command to adjust the hydraulic throttling element of the heave compensation device and gradually locking the heave compensation hydraulic cylinder of the heave compensation device to complete the fault self-protection includes the following steps: If the instability prediction index exceeds the preset risk threshold, the target damping coefficient is calculated using a fuzzy logic controller based on the extent to which the instability prediction index exceeds the risk threshold. The target damping coefficient is converted into an opening control signal for the hydraulic throttling element in the heave compensation device, and the proportional throttling valve of the heave compensation device is driven to reduce the flow area. During the adjustment process, the movement deceleration of the hydraulic cylinder of the heave compensation device is monitored in real time. If the movement deceleration exceeds the safety limit, the current throttling opening of the proportional throttle valve is maintained until the hydraulic cylinder completely stops moving and locks.

[0013] In a second aspect, the present invention also provides a fault prediction and self-protection system for diving bell heave compensation, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the fault prediction and self-protection method for diving bell heave compensation as described in any one of the first aspects.

[0014] Thirdly, the present invention also provides a computer-readable storage medium storing instructions, characterized in that, when executed by a processor, the instructions cause the processor to be configured to perform a fault prediction and self-protection method for diving bell heave compensation according to any one of the first aspects.

[0015] The beneficial effects of this invention are: This invention maps the relationship between hydraulic cylinder displacement and accumulator pressure into a hysteresis feature field image and performs geometric moment analysis, enabling quantitative characterization of the deterioration trend of the dynamic response characteristics of the heave compensation device at the system-wide level. By performing homomorphic deconvolution and cepstral filtering on the pressure signal, friction texture signals that directly reflect the mechanical contact state of internal seals and guides can be separated from strong background noise. By calculating their energy entropy values, accurate quantification of the deterioration of micro-friction states is achieved. Next, a prediction framework based on a nonlinear cusp catastrophe model is constructed, mapping macroscopic topological feature descriptors to state points in the control variable space, and calculating the geometric distance from this point to the boundary of the model bifurcation set. The set distance provides a quantifiable measure of the critical point at which the system is close to sudden failure. Finally, by combining this geometric distance with the micro-energy entropy value, a highly reliable instability prediction index is calculated, upgrading fault diagnosis to the level of early prediction and active self-protection against sudden system instability failure. This solves the technical problem that traditional methods struggle to provide early warning of sudden system instability failure caused by multi-factor coupling. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the fault prediction and self-protection method for diving bell heave compensation in one embodiment of this application.

[0017] Figure 2 This is a schematic diagram of a binary image of a hysteresis feature field image in one embodiment of this application.

[0018] Figure 3 This is a schematic diagram of a cepstral energy probability distribution histogram in one embodiment of this application. Detailed Implementation

[0019] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0020] The terms "first," "second," etc., used in the specification and claims of this application are used to distinguish similar objects and not to describe a specific order or sequence. It should be understood that such use of data can be interchanged where appropriate so that embodiments of this application can be implemented in orders other than those illustrated or described herein, and the objects distinguished by "first," "second," etc., are generally of the same class and the number of objects is not limited; for example, a first object can be one or more. Furthermore, in the specification and claims, "and / or" indicates at least one of the connected objects, and the character " / " generally indicates that the preceding and following objects are in an "or" relationship.

[0021] Figure 1 This is a flowchart illustrating a fault prediction and self-protection method for diving bell heave compensation in one embodiment. It should be understood that, although... Figure 1 The steps in the flowchart are shown sequentially as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified herein, there is no strict order in which these steps are executed, and they can be performed in other orders. Figure 1 At least some steps in the process may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps. For example Figure 1 As shown, the fault prediction and self-protection method for diving bell heave compensation disclosed in this invention specifically includes the following steps: S101. Obtain the real-time operating data of the heave compensation device, extract the hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle from the real-time operating data, map them to a two-dimensional coordinate system to construct a closed hysteresis loop, and convert the closed hysteresis loop into a binary hysteresis feature field image.

[0022] The system synchronously acquires displacement signals from the hydraulic cylinders and pressure signals from the accumulator cavity in the heave compensation device using high-frequency sensors. These two physical quantities represent the system's motion state and dynamic response, respectively, and their coupling accurately reflects the nonlinear friction and energy dissipation within the mechanical system. Since the original discrete data points may exhibit slight asynchrony or sampling interval fluctuations on the time axis, a cubic spline interpolation algorithm is first used to reconstruct the displacement and pressure data within a single heave cycle, generating a smooth and continuous closed trajectory, effectively avoiding feature extraction errors caused by data breakpoints. Subsequently, the processed displacement and pressure data are mapped to the horizontal and vertical axes of a two-dimensional Cartesian coordinate system, constructing a closed hysteresis loop that intuitively characterizes the system's damping properties. To convert this analog signal into a digital image format that can be processed by computer vision algorithms, a two-dimensional grid resolution adapted to the system's range is set, and the coordinate values ​​of the hysteresis loop are normalized to a unit interval and projected onto this grid. Grid positions traversed by the trajectory are assigned foreground pixels, and other positions are assigned background pixels, thus generating a binarized initial image. Considering that signal noise may cause minor breaks in the trajectory, morphological closing operations are applied to fill holes and smooth the image, ultimately obtaining a hysteresis feature field image that can serve as a benchmark for topological feature extraction. This completes the dimensionality reduction and feature solidification process from one-dimensional time series to two-dimensional topological graphics.

[0023] S102. Perform geometric moment analysis on the hysteresis feature field image and calculate the topological feature descriptors for the quantized hysteresis loop discreteness and skewness features.

[0024] The process involves calculating the zero-order moment of the image to obtain the closed area enclosed by the hysteresis loop, which physically corresponds to the total energy dissipated by the system due to friction and damping in a single cycle. Simultaneously, the first-order moment is calculated to determine the geometric centroid position of the image in the coordinate system, reflecting the equilibrium center of the system during reciprocating motion. To eliminate the interference of absolute positional changes in the image within the coordinate system on feature analysis, second- and third-order central moments with smoothness and translation invariance are calculated using the geometric centroid. Based on this, specific invariant feature formulas are constructed by normalizing the central moments. The combination of second-order central moments is defined as the discreteness feature, which physically describes the inertial distribution of the hysteresis loop along the principal axis, i.e., the thickness of the image, directly related to the clearance changes caused by wear of the hydraulic cylinder seals. The combination of third-order central moments is defined as the skewness feature, which quantifies the asymmetry of the image relative to the centroid and can detect uneven wear or abnormal lateral forces occurring on the piston rod.

[0025] S103. Homomorphic deconvolution is used to process the accumulator pressure data and transform the accumulator pressure data to the cepstrum domain. Low-frequency components of the accumulator pressure data are filtered out by inverted high-pass filtering. At the same time, the friction texture signal characterizing the mechanical contact state is separated and the energy entropy value of the friction texture signal is calculated.

[0026] Since the pressure signal in the time domain is a convolution of the system's inherent characteristics and external excitation, direct filtering is difficult to separate it. Therefore, a Hamming window is first applied to the data to reduce spectral leakage. A Fast Fourier Transform (FFT) is then used to transform it to the frequency domain, and the natural logarithm of the spectral amplitude is taken. This nonlinear operation transforms the original multiplicative convolution relationship into an additive one. Subsequently, an Inverse Fast Fourier Transform (IFFT) is performed to map the signal to the inverse frequency domain. In the inverse frequency domain, the slowly varying wave component is located in the low inverse frequency range, while the rapidly varying mechanical friction characteristics are located in the high inverse frequency range. An optimal inverse frequency cutoff threshold is set according to the sea state cycle, and an inverse high-pass filter is designed to filter out low-frequency interference, retaining only the high inverse frequency residual component. An inverse homomorphic transform is then performed on this residual component to restore it to the time domain, allowing the extraction of the pure friction texture signal. Finally, the signal is treated as a random process to calculate the probability density function, and the energy entropy value is calculated using the Shannon entropy formula. The entropy value directly quantifies the disorder and complexity of the mechanical contact interface; a higher entropy value indicates more severe wear or erosion of the friction surface.

[0027] S104. Map the discreteness characteristics to splitting factor control variables and the skewness characteristics to normal factor control variables, construct a control variable space composed of splitting factor control variables and normal factor control variables, and determine the current state point of the heave compensation device in the control variable space.

[0028] The acquisition system collects a large amount of historical data under healthy baseline conditions to establish statistical distribution models for discreteness and skewness characteristics. The discreteness and skewness characteristics calculated at the current moment are substituted into these models to calculate standardized Z-scores, eliminating dimensional differences and determining the degree of deviation from normal values. Then, using the Sigmoid and hyperbolic tangent functions as compression mapping tools, these two unbounded statistical scores are transformed into bounded control variables: the discrete characteristics are mapped to a splitting factor, which controls the opening size of the system's potential energy surface, i.e., the potential probability of a sudden change in the system; the skewness characteristics are mapped to a normality factor, which determines the tilt direction of the system on the potential energy surface, i.e., the direct driving force triggering a sudden change. This constructs a two-dimensional control variable space, where each operating state of the heave compensation device is positioned as a specific current state point.

[0029] S105. Calculate the geometric distance from the current state point to the boundary of the bifurcation set of the nonlinear cusp catastrophe model based on the coordinates of the current state point in the control variable space, and calculate the instability prediction index by combining the geometric distance and the energy entropy value.

[0030] According to the geometric definition of cusp catastrophe theory, there exists a bifurcation set in the control variable space, which is the critical boundary for the system's potential energy function to undergo a topological abrupt change. Mathematically, it is represented by a cusp-shaped curve. Crossing this curve means that the system will transition from a stable monostable state to an uncontrollable bistable state or experience violent oscillations. An iterative search algorithm is used to accurately calculate the shortest Euclidean geometric distance from the current state point to the boundary of the bifurcation set curve. This distance directly reflects the system's stability margin; the smaller the distance, the greater the risk of abrupt failure. Since the geometric distance alone does not include microscopic contact state information, the friction texture energy entropy value is introduced as a correction term. After normalizing the energy entropy value, a weighting coefficient is set and it is fused with the geometric distance to calculate the final instability prediction index. The instability prediction index considers not only the topological stability of the macroscopic motion trajectory but also the degree of degradation of the microscopic friction interface, thus achieving high-precision prediction of the failure trend of the heave compensation device.

[0031] S106. When the instability prediction index meets the preset risk threshold condition, a variable damping control command is generated to adjust the hydraulic throttling element of the heave compensation device, and the heave compensation hydraulic cylinder of the heave compensation device is gradually locked to complete the fault self-protection.

[0032] When the calculated instability prediction index exceeds the preset safety risk threshold, it indicates that the heave compensation device is on the verge of functional failure or mechanical damage, immediately triggering a graded self-protection control strategy. Instead of abruptly applying emergency braking to prevent pipeline rupture or secondary injury to personnel inside the diving bell due to hydraulic shock, the control core generates variable damping control commands. Specifically, a fuzzy logic controller calculates the target damping coefficient required for the current operating condition based on the magnitude of the instability index exceeding the threshold, and converts this coefficient into a current signal to drive the proportional throttle valve. By rapidly and smoothly adjusting the valve opening of the hydraulic throttle element, the return oil flow area is gradually reduced, thereby applying a gradually increasing hydraulic damping force to the heave compensation hydraulic cylinder. During this process, the system monitors the cylinder's speed and deceleration in real-time in a closed loop, ensuring that the deceleration remains within permissible safety limits. As the fluid resistance increases smoothly, the piston movement of the hydraulic cylinder is gently suppressed and eventually decelerates to a complete stop. Once the hydraulic cylinder speed is detected to be zero, the system immediately locks the mechanical locking mechanism or closes the shut-off valve, safely locking the device in the current position, completing the self-protection process from fault warning to active flexible shutdown.

[0033] In one embodiment, extracting hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle from real-time operational data, mapping them to a two-dimensional coordinate system to construct a closed hysteresis loop, and converting the closed hysteresis loop into a binary hysteresis feature field image includes the following steps: Set a two-dimensional grid resolution that is adapted to the displacement range of the hydraulic cylinder and the pressure variation range of the accumulator, and initialize the two-dimensional grid. The two-dimensional grid resolution is used to define the pixel density of the hysteresis feature field image. The hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle in the real-time operation data are processed by cubic spline interpolation to generate a smooth and continuous closed hysteresis loop trajectory. The coordinates of the closed hysteresis loop trajectory are normalized to a unit interval, and the normalized coordinates are mapped to a two-dimensional grid. The grid cells passed by the closed hysteresis loop trajectory are marked as foreground pixels, and the grid cells not passed by the closed hysteresis loop trajectory are marked as background pixels, thus generating an initial two-dimensional grid image. The initial two-dimensional mesh image is preprocessed using morphological closing operations to obtain a two-dimensional mesh image. The two-dimensional mesh image is then transformed into a binary matrix containing only 0 and 1 values ​​to generate a hysteresis feature field image for feature extraction.

[0034] In this embodiment, refer to Figure 2 Based on the maximum stroke range of the hydraulic cylinder in the diving bell heave compensation device and the highest possible pressure peak of the accumulator, the physical boundaries of the two-dimensional feature field image are determined. A high-precision two-dimensional grid resolution is set. ,in The number of mesh divisions representing the lateral displacement axis. The number of grid divisions representing the longitudinal pressure axis directly defines the pixel density and detail resolution of the final generated hysteresis feature field image. A zero-matrix is ​​initialized as the underlying digital canvas; the dimensions of this matrix strictly correspond to the set grid resolution, and each element represents a tiny rectangular region in the coordinate system. Within a single heave cycle, cubic spline interpolation is applied to reconstruct the collected hydraulic cylinder displacement and accumulator pressure sequences. This algorithm uses piecewise cubic polynomials to construct smooth curves between every two adjacent data points, ensuring the continuity of the first and second derivatives at connection points, thus mathematically eliminating the jagged effect caused by the sampling interval. By inserting high-density virtual sampling points between the original time steps, smoother and more continuous displacement and pressure interpolation sequences are generated. These two sequences correspond to each other on the two-dimensional plane, jointly depicting a closed hysteresis loop trajectory, which completely records the energy dissipation path of the system from compression to rebound.

[0035] After obtaining a smooth, continuous trajectory, the interpolated displacement and pressure data are normalized to linearly compress all values ​​to a unit interval of [0,1], eliminating the interference of different physical dimensions on image generation. Then, the normalized coordinate values ​​are multiplied by the grid resolution and rounded down to calculate the specific index position of each point on the trajectory in the two-dimensional grid matrix. All sampling points on the trajectory are traversed, and the element value at the corresponding index position in the matrix is ​​changed from the initial 0 to 1, marked as a foreground pixel, representing that the grid cell is traversed or occupied by the physical trajectory; while grid cells not touched by the trajectory remain at 0, marked as background pixels. Due to sensor noise or extremely small discontinuities during the interpolation process, the initially generated grid image may contain minute breaks or isolated noise holes within the trajectory lines, which can interfere with subsequent connected component analysis and topological feature extraction. Therefore, it is necessary to introduce the closing operation processing technique from digital image morphology. Specifically, firstly, a structuring element of a preset size is slid across the image matrix to perform a dilation operation on the edges of foreground pixels, filling small cracks and connecting adjacent breakpoints. Then, the same structuring element is used to perform an erosion operation, restoring the image contour to its original thickness while preserving the connection effect and removing false edges introduced by the dilation. Finally, it is ensured that the image matrix is ​​strictly converted into a binary matrix containing only logic 0 (background) and logic 1 (foreground).

[0036] In one embodiment, performing geometric moment analysis on the hysteresis feature field image and calculating the topological feature descriptor for quantizing the discreteness and skewness features of the hysteresis loop includes the following steps: The zeroth moment of the hysteresis feature field image is calculated to obtain the closed area of ​​the hysteresis loop, and the geometric centroid coordinates of the hysteresis loop in the two-dimensional coordinate system are determined by calculating the first moment of the hysteresis feature field image. The second and third central moments of the hysteresis feature field image are calculated based on the geometric centroid coordinates. The second and third central moments are then normalized using the zero-order moment to obtain the normalized central moments. The first invariant characteristic formula is constructed using the second-order central moment in the normalized central moments, and the sum of the moments of inertia in the orthogonal directions is defined as the discreteness characteristic that characterizes the discreteness of the figure. The second invariant feature formula is constructed using the third central moment in the normalized central moments. The sum of the projection skewness in the orthogonal directions is defined as the skewness feature that characterizes the degree of skewness of the figure. The calculated discreteness feature and the skewness feature are combined to form a topological feature descriptor.

[0037] In this embodiment, the binarized hysteresis feature field image is treated as a two-dimensional discrete density distribution function, and its basic geometric properties are extracted using the principle of image moments. Image moments are weighted statistical features of image pixel intensity and position coordinates, where the zeroth moment represents the total quality of the image, equivalent to the sum of foreground pixels in a binary image. By traversing the entire image matrix and accumulating the gray values ​​of all pixels, the zeroth moment of the hysteresis feature field image is calculated. This value directly corresponds to the geometric area enclosed by the closed hysteresis loop at the physical level, reflecting the total energy dissipated by the hydraulic system due to damping and friction during a single heave cycle. Based on obtaining the zero-order moment, coordinate position information is further introduced to calculate the first-order moment of the image. The weighted sum of pixel coordinate values ​​in the horizontal and vertical directions is calculated separately to obtain the horizontal first-order moment. With the first perpendicular moment These two parameters quantify the moment distribution of the image about the coordinate axes. This is achieved using the ratio of the first moment to the zeroth moment, i.e. and It can calculate the geometric centroid coordinates of the hysteresis loop in a two-dimensional pixel coordinate system. The geometric centroid coordinates do not represent the average operating point or equilibrium position of the heave compensation device under current operating conditions. The coordinates of each foreground pixel in the image are centered based on the geometric centroid coordinates, i.e., the original coordinates are subtracted from the centroid coordinates. The second and third central moments of the image are then calculated using the transformed relative coordinates. Second central moment... , and The graph's distribution dispersion about the horizontal axis, vertical axis, and diagonal axis, respectively, and the third-order central moment are described. , These capture the asymmetry and degree of distortion of the graphic shape.

[0038] Because changes in image size or resolution cause overall scaling of the central moments, normalization is necessary to eliminate the interference of this scale factor. This is achieved using the zeroth moment. As a normalization factor, according to the formula By performing division operations on the central moments of each order, the normalized central moments are obtained. Here, p and q represent the order of moments in the x and y directions, respectively. The obtained normalized central moments not only eliminate the effects of image translation but also the effects of image scaling, extracting the essential features that purely reflect the topological shape and structure of the hysteresis loop. In image analysis, the second-order moments mainly reflect the mass distribution of an object around its center of mass. By combining the second-order normalized central moments in orthogonal directions, a first invariant feature formula is constructed to calculate the discreteness feature. The discreteness characteristic is physically equivalent to the polar rotational inertia of an image about its centroid, describing the average degree of diffusion of image pixels relative to the geometric center.

[0039] In the fault diagnosis of heave compensation devices, the morphological changes of the hysteresis loop are directly related to the system's damping characteristics and leakage state. When the hydraulic cylinder's sealing performance deteriorates or friction increases, the hysteresis loop typically becomes wider or more divergent, causing the image pixels to be distributed further away from the centroid in space, thus significantly increasing the calculated discreteness feature value. Conversely, if the system operates compactly and efficiently, the loop is relatively narrow and the discreteness feature value is smaller. Therefore, the discreteness feature, as a highly condensed scalar indicator, can effectively characterize the discreteness of the hysteresis loop, intuitively reflecting the severity of energy dissipation and potential mechanical clearance problems in the system, thus transforming the complex two-dimensional graphical form into a single measurable fault monitoring variable.

[0040] In statistics, the third moment corresponds to skewness, which can keenly reflect the direction and degree of the tailing of a distribution. Using the third component of the normalized central moment, a second invariant characteristic formula is constructed to calculate the degree of skewness. Specifically, by combining the projected skewness in the horizontal and vertical directions, the formula is used... This formula defines the degree of skewness in the graph. It integrates the third-order projection components of the image along the orthogonal axes, quantifying the degree of distortion or tilt of the hysteresis loop relative to the central axis. In actual operation, if the heave compensation hydraulic cylinder experiences lateral wear, the guide mechanism jams, or the sensor zero point drifts, the hysteresis loop will often exhibit obvious asymmetry or local bulges. The magnitude of the numerical value directly corresponds to the severity of this asymmetric distortion. Finally, the calculated discreteness features and skewness features are vectorized and combined to form a topological feature descriptor containing two dimensions. The topological feature descriptor fully encodes the macroscopic contour and microscopic distortion information of the hysteresis loop.

[0041] In one embodiment, the process of processing the accumulator pressure data using homomorphic deconvolution and transforming the accumulator pressure data to the cepstrum domain, filtering out the low-frequency components of the accumulator pressure data using an inverted high-pass filter, and simultaneously separating the friction texture signal characterizing the mechanical contact state and calculating the energy entropy value of the friction texture signal includes the following steps: A Hamming window is applied to the collected accumulator pressure data, and a fast Fourier transform is performed on the windowed accumulator pressure data to obtain the pressure spectrum. Taking the natural logarithm of the amplitude of the pressure spectrum transforms the convolutional relationship in the time domain of the pressure spectrum into an additive relationship in the logarithmic amplitude spectrum in the frequency domain. Perform an inverse fast Fourier transform on the logarithmic amplitude spectrum in the frequency domain to map the accumulator pressure data to the inverted frequency domain to obtain the cepstrum sequence of the accumulator pressure data; The inverted frequency cutoff threshold is set according to the wave compensation period of the heave compensation device, and an inverted high-pass filter is constructed to filter out the low inverted frequency components representing wave motion and fluid pressure in the cepstral sequence. An inverse homomorphic transform is performed on the high inverse frequency residual component after processing by the inverse high-pass filter, and the high inverse frequency residual component is restored to the time domain to obtain the friction texture signal; The probability density function of the friction texture signal is calculated, and the entropy value of the probability density function is calculated based on the Shannon entropy formula to obtain the energy entropy value that quantifies the complexity of the mechanical contact surface.

[0042] In this embodiment, before performing time-domain signal conversion, a Hamming window function w(n) is first applied to the acquired raw discrete data sequence of accumulator pressure x(n). The Hamming window can effectively reduce sidelobe amplitude and reduce spectral interference. By multiplying the raw signal with the window function point by point, the windowed signal sequence is obtained. The windowed signal was then processed using the Fast Fourier Transform algorithm, mapping it from the time domain to the frequency domain, and the complex form of the pressure spectrum was calculated. The pressure spectrum contains the amplitude and phase information of the signal at various frequency points. In the time domain, the pressure signal of the heave compensation device is actually a complex signal formed by the convolution of external excitation (such as wave motion) and the system's own transmission characteristics (such as mechanical friction and pipeline damping). This convolution relationship makes it difficult for simple linear filtering to completely separate the two; however, this nonlinear convolution effect can be decoupled using a logarithmic transformation step.

[0043] Specifically, first, the amplitude spectrum |X(k)| of the pressure spectrum is calculated, and then the natural logarithm of this amplitude spectrum is taken to obtain the logarithmic amplitude spectrum. This transforms the multiplicative or convolutional combination relationship of the original signal in the time domain into a linear additive relationship in the logarithmic spectrum of the frequency domain. After the transformation, the originally entangled low-frequency wave signal components and high-frequency mechanical friction signal components become simple superposition terms in the logarithmic spectrum. Next, the inverse fast Fourier transform algorithm is used to process the logarithmic amplitude spectrum to calculate the cepstrum sequence of the accumulator pressure data. In the reciprocal frequency domain, the horizontal axis represents the reciprocal frequency with the dimension of time. Different signal components are mapped to different positions on the reciprocal frequency axis according to the rate of change of their spectra: components with slow spectral changes correspond to the envelope of the signal or low-frequency excitation and are concentrated in the low reciprocal frequency region, while components with rapid spectral changes correspond to the fine structure or high-frequency texture of the signal and are distributed in the high reciprocal frequency region.

[0044] Considering that the primary operating context of heave compensation devices is to cope with ocean wave fluctuations, their wave compensation period is typically between a few seconds and tens of seconds, which corresponds to a relatively low numerical range in the cepstral domain. Based on actual marine environmental parameters and the device's design specifications, a clear cepstral cutoff threshold is determined. This threshold distinguishes between low-frequency wave interference components and high-frequency mechanical characteristic components. A high-pass filter is constructed based on this threshold, with its transfer function designed as either a step function or a smooth transition function. All sequence values ​​with cepstral frequencies less than the threshold are zeroed or significantly attenuated, thus filtering out low-frequency components representing wave motion and macroscopic fluid pressure. The retained cepstral sequence mainly consists of high-frequency components, which physically correspond to rapidly changing details in the pressure signal, i.e., signal fluctuations caused by mechanical friction and changes in surface roughness.

[0045] After high-pass filtering, a residual cepstrum sequence containing only high-frequency components is obtained. A Fast Fourier Transform (FFT) is performed on this filtered high-frequency cepstrum sequence to return it to the frequency domain, yielding the high-frequency components of the logarithmic amplitude spectrum. Next, an exponential function is used to remove the influence of the logarithmic operation, restoring the linear amplitude spectrum, and the complex spectrum is reconstructed by combining it with the phase information of the original signal. Finally, an inverse FFT is performed on the reconstructed spectrum to map the signal from the frequency domain back to the time domain, obtaining the friction texture signal. This restored time-domain signal removes large-amplitude wave pressure fluctuations, retaining only the high-frequency vibration texture related to mechanical friction and the contact state of the seals. This texture reflects the contact condition between the hydraulic cylinder piston rod and the seals, such as whether there is dry friction, scratches, or lubrication failure. Finally, the extracted time-domain friction texture signal is treated as a random process, and its amplitude distribution is statistically analyzed. The signal amplitude range is divided into several discrete intervals, and the number of signal data points falling within each interval is counted, thus calculating the probability density function of the friction texture signal, denoted as p(x). Based on this probability distribution, the energy entropy value is calculated using the Shannon entropy formula. Energy entropy represents the uncertainty or richness of information in a signal. In heave compensation devices, when the mechanical surface is smooth and well-lubricated, the friction signal is regular and uniform, with a low entropy value. However, when wear, erosion, or lubrication failure occurs, the friction signal becomes chaotic, containing more random mutations, leading to a significant increase in entropy. Therefore, the calculated energy entropy value directly quantifies the complexity and disorder of the mechanical contact surface.

[0046] In one embodiment, performing an inverse homomorphic transform on the high-frequency inverse residual component after processing with the inverse high-pass filter to restore the high-frequency inverse residual component to the time domain to obtain the friction texture signal includes the following steps: An inverse complex cepstral transform is performed on the high-frequency inverse residual component after processing by the inverse high-pass filter to restore the high-frequency inverse residual component to the time domain, generating a mixed time-domain sequence containing mechanical friction and impact characteristics and background Gaussian white noise. An impulse response filter with adjustable filter coefficients is constructed. The kurtosis value of the output signal of the impulse response filter is defined as the objective optimization function. The kurtosis value is used to characterize the significance of non-Gaussian impulse components in the signal. Stick-slip vibration caused by mechanical friction has high kurtosis and non-Gaussian properties. The minimum entropy deconvolution algorithm is used to iteratively update the adjustable filter coefficients of the impulse response filter until the objective optimization function converges to the maximum value, thereby generating an inverse filter to maximize the enhancement of friction and impact characteristics. The convolution of the mixed time-domain sequence is performed using the converged inverse filter to eliminate the signal diffusion effect caused by the transmission path and suppress the background Gaussian white noise, thereby reconstructing the friction texture signal.

[0047] In this embodiment, for the retained high-frequency reciprocal residual components, an inverse complex cepstral transform technique is first applied to perform a forward Fourier transform on the residual components to return them to the logarithmic frequency domain. Then, exponential operations are used to eliminate the influence of the previous logarithmic processing, thereby restoring the linearity of the amplitude spectrum. Finally, an inverse Fourier transform is performed in conjunction with the original phase information to map the signal back to the time domain. A finite impulse response filter is then constructed, which has a set of dynamically adjustable coefficient vectors f of length L. The filter reshapes the input signal by adjusting its frequency response characteristics. To guide the adjustment direction of the filter, an objective optimization function is defined, specifically using the kurtosis value of the filter's output signal y. Kurtosis is a dimensionless fourth-order statistic, and its calculation formula is... ,in The mean, The standard deviation is denoted by . This index is exceptionally sensitive to impulsive components in the signal: for pure Gaussian noise, its kurtosis value is close to 3; while for signals containing significant impulsive or pulse components, its kurtosis value will be significantly greater than 3. Because the stick-slip phenomenon caused by seal wear or poor lubrication in the heave compensation device generates sharp non-Gaussian transient vibrations, these vibrations have extremely high kurtosis characteristics. Therefore, setting the kurtosis value of the output signal as the optimization objective is essentially finding a filtering method that ensures the filtered signal contains as many and as clearly as possible the mechanical impulsive components.

[0048] After establishing the optimization objective, the minimum entropy deconvolution algorithm is used to specifically solve for the optimal filter coefficients. Specifically, an inverse filter can be designed to counteract the tailing or diffusion effect caused by the system transmission path on the original impulse signal, thus refocusing the originally dispersed impulse energy. The algorithm initializes a unit impulse as the initial coefficients of the filter and then enters a loop iteration: in each iteration, the update gradient of the filter coefficients is calculated based on the higher-order moment characteristics of the current output signal. Matrix operations are used to continuously correct the filter coefficient vector, causing the kurtosis value of the output signal to gradually increase. When the change in kurtosis value calculated in two consecutive iterations is less than a preset small threshold, or when the maximum number of iterations is reached, the objective function is considered to have converged to a global or local maximum. At this point, the filter coefficients are locked, generating the final inverse filter. Next, an inverse filter is used to perform convolution processing on the mixed time-domain sequence. At the physical level, this process produces two significant effects: First, the inverse filter effectively compensates for waveform distortion caused by structural damping and multipath effects during signal propagation from the friction source to the sensor, recompressing the originally stretched and blurred impact waveform into sharp, clear pulses and eliminating signal diffusion. Second, because the filter is designed to maximize non-Gaussianity, it naturally suppresses background white noise that follows a Gaussian distribution, as the noise has low kurtosis and cannot pass the filter's filtering. Therefore, the signal reconstructed by the inverse filter's convolution processing is a high-fidelity friction texture signal, in which the impact characteristics caused by mechanical friction are significantly enhanced, and the signal-to-noise ratio is significantly improved.

[0049] In one embodiment, setting a cepstral frequency cutoff threshold based on the wave compensation period of the heave compensation device and constructing a high-pass cepstral filter to filter out low-frequency components representing wave motion and fluid pressure in the cepstral sequence includes the following steps: Calculate the amplitude envelope of the cepstrum sequence and normalize the amplitude envelope to construct a histogram of the cepstrum energy probability distribution; The information complexity of the distribution histogram is quantified by Raney entropy, and the sum of the Raney entropy of the background class cepstral distribution and the foreground class cepstral distribution is defined as the objective function of the candidate cepstral cutoff threshold. Traverse all candidate cepstral frequency cutoff thresholds within the domain of the cepstral sequence, solve the global maximum value of the objective function using a genetic algorithm, and determine the candidate cepstral frequency cutoff threshold corresponding to the global maximum value as the optimal cepstral frequency cutoff threshold; Based on the optimal cepstral frequency cutoff threshold, a cosine tapered window function with smooth transition characteristics is constructed. The cosine tapered window function is multiplied with the cepstral sequence to suppress wave fluid interference at low cepstral frequencies while preserving the friction characteristics of high cepstral frequencies.

[0050] In this embodiment, the absolute value of the cepstral sequence is calculated to obtain the amplitude envelope |c(n)|, which reflects the energy strength of different cepstral components. Since the original amplitudes may span large orders of magnitude, direct statistical analysis would face numerical instability; therefore, the amplitude envelope is normalized. This is achieved by dividing the amplitude at each point by the sum of the amplitudes of the entire sequence. Converting physical energy values ​​into mathematical probability values Based on these probability values, a distribution histogram reflecting the characteristics of the cepstral energy distribution is constructed, such as... Figure 3 As shown in the histogram, the horizontal axis represents the index position of the inverted frequency, and the vertical axis represents the energy probability density corresponding to that position.

[0051] After obtaining the energy distribution histogram, for any given candidate cepstral cutoff threshold t, the entire cepstral distribution is divided into two parts: one part is the background distribution representing low-frequency disturbances such as wave motion, and the other part is the foreground distribution representing high-frequency features such as mechanical friction. The normalized probabilities of these two distributions are calculated separately, and then the formula is used... Calculate their respective Reni entropy, where It is the sum of probabilities for corresponding categories. Define an objective function. The objective function represents the sum of the Raney entropy of the background class and the foreground class. According to the principle of maximum entropy, when the system is in the most reasonable segmentation state, the sum of the information content of the two parts should reach its maximum, that is, the uniformity within the two distributions is the highest, and the difference between the categories is the greatest.

[0052] Next, a genetic algorithm is used to solve for the globally optimal threshold. Specifically, an initial candidate threshold population is randomly generated within the domain of the cepstrum, with each threshold encoded as a chromosome. Then, an iterative evolutionary process begins: the objective function value corresponding to each individual is calculated as its fitness; a selection mechanism is used to weed out the less fit individuals, retaining those with high fitness; crossover and mutation operations are performed on the selected individuals to generate a new offspring population. After multiple generations of evolution, the population gradually converges towards the peak region of the objective function. Finally, when the evolution reaches a preset number of generations or the fitness no longer significantly improves, the individual with the highest fitness in the current population is selected and decoded as the optimal cepstrum cutoff threshold.

[0053] After determining the optimal cepstral cutoff threshold using a genetic algorithm, a cosine taper window function (such as a Tukey window or a similar smooth high-pass window) with smooth transition characteristics is constructed. The cosine taper window function takes a value of 0 or rapidly decays to 0 in the low cepstral frequency range from 0 to the optimal cepstral cutoff threshold to completely suppress wave components; in the transition region after the optimal cepstral cutoff threshold, it smoothly rises to 1 according to a cosine curve; and in the high cepstral frequency range, it remains at 1 to fully pass and preserve friction characteristics. The constructed window function sequence is then multiplied point-by-point with the original cepstral sequence c(n). This soft-thresholding method thoroughly filters out high-energy low-frequency interference representing fluid pressure in the frequency domain while maximally preserving the high-frequency weak signals representing mechanical states. Simultaneously, the smooth transition edges effectively avoid time-domain ringing during signal reconstruction.

[0054] In one implementation, mapping the discreteness characteristics to splitting factor control variables and the skewness characteristics to normality factor control variables, constructing a control variable space composed of splitting factor control variables and normality factor control variables, and determining the current state point of the heave compensation device in the control variable space includes the following steps: A nonlinear cusp mutation model of the heave compensation device is constructed and a standard cusp mutation potential function is established. The standard cusp mutation potential function includes a state variable, a quartic term, a quadratic term weighted by a splitting factor, and a linear term weighted by a normality factor. Historical operating data of the heave compensation device under healthy conditions were collected, and a baseline distribution model was established by calculating the dispersion and skewness characteristics of the historical operating data. Calculate the first standardized Z-score of the discreteness feature extracted at the current time relative to the baseline distribution model, and the second standardized Z-score of the skewness feature extracted at the current time relative to the baseline distribution model; The first standardized Z-score is compressed and mapped using a nonlinear Sigmoid function to generate a splitting factor control variable corresponding to the quadratic term of the standard cusp mutation potential function. The hyperbolic tangent function is used to compress and map the second standardized Z-score, generating a normal factor control variable corresponding to the first term of the standard cusp mutation potential energy function; A control variable space consisting of splitting factor control variables and normal factor control variables is constructed, and the current state point of the heave compensation device in the control variable space is determined based on the standard cusp mutation potential energy function.

[0055] In this embodiment, the nonlinear cusp catastrophe model can qualitatively capture the state transition behavior of complex systems using polynomials. The total potential energy function V(x) of the heave compensation device control system is defined, where x is a single state variable representing the system behavior. The standard cusp catastrophe potential energy function is constructed as follows: In this equation, the quartic term... The fundamental constraints of the potential energy trap are provided, ensuring the global stability of the system; quadratic terms The potential energy curve is weighted by the parameter u and called the splitting factor. It determines whether the potential energy curve has a single valley or a double valley, that is, whether the system has the possibility of multiple stable states. The first-order term vx is weighted by the parameter v and called the normality factor. The normality factor causes the potential energy curve to tilt to the left and right, driving the asymmetric shift of the system state.

[0056] Next, historical sensor data marked as healthy or operating normally needs to be selected from the heave compensation device's maintenance database. This data is then processed using the aforementioned method to extract the dispersion and skewness characteristics for each moment in batches. Statistical methods are then used to analyze the large-sample distribution patterns of these two characteristics, and the mean value of the dispersion characteristic under healthy conditions is calculated. with standard deviation and the mean of the skewness characteristics. with standard deviation These two sets of statistical parameters constitute the system's baseline distribution model. During real-time monitoring, standardization is necessary to eliminate the significant differences in units and numerical ranges among different physical characteristics. This standardization process is crucial when extracting new discrete characteristics from the current operational data. and skewness characteristics Then, Z-score standardization is performed using the parameters of the established baseline distribution model. The specific calculation formula is: First Standardized Z-score and the second standardized Z-score These two Z-scores represent how many standard deviations the current feature value deviates from the healthy mean. This reflects the significance of the changes in the plumpness or thinness of the hysteresis loop. This reflects the significant asymmetry and distortion of the graphic.

[0057] Since the control variables u and v in cusp catastrophe models are typically defined within specific bounded intervals, while the statistically obtained Z-score is theoretically unbounded, a nonlinear compression mapping is required. For the first standardized Z-score, representing the degree of discreteness, a sigmoid function is used for transformation to generate the splitting factor control variable u. The mapping formula is designed as follows: ,in The adjustment coefficient, u, is used to set the range of values ​​for u. The sigmoid function's S-shaped curve characteristic makes the mapping process linear and gradual when the Z-score is small, while it tends to saturate at extremely large values, effectively preventing the model parameters from overflowing due to drastic fluctuations in discrete characteristics. Physically, as the dispersion of the hysteresis loop increases, the value of u changes, driving the system's potential energy surface to gradually split from a single stable point into two potential stable points, simulating the multistable or unstable trend of the system caused by the increase of wear gaps.

[0058] Similarly, for the second standardized Z-score representing the degree of skewness, a hyperbolic tangent function is used for mapping to generate a normality factor control variable v. The calculation formula is as follows: ,in The adjustment coefficient is used. The hyperbolic tangent function smoothly maps the input value to... Within the interval, it exhibits good linearity near the origin. In the catastrophe model, v represents the external driving force that disrupts the system's symmetry. When the hysteresis loop shows significant skewness distortion, the mapped v value deviates significantly from zero, causing the potential energy curve to tilt, making the system more prone to instability. The introduction of the hyperbolic tangent function not only normalizes the numerical range but also simulates the nonlinear response characteristics of the system to skewness disturbances. Small skewnesses may be absorbed by the system's own stiffness, while large skewnesses quickly transform into driving forces for catastrophes, specifically describing the destructive effect of asymmetric faults on system stability.

[0059] After the above mapping, the splitting factor u and the normality factor v are obtained. These two variables together constitute a two-dimensional Cartesian coordinate system, i.e., the control variable space. The operating state of the heave compensation device at any given time can be uniquely located as a coordinate point P(u,v) in this space, which is the current state point. Combined with the standard cusp catastrophe potential energy function, this state point is not only a geometric location but also contains dynamic meaning. On the (u,v) plane, different regions correspond to different topological structures of the potential energy function: some regions have only one minimum value (stable), while others have two minimum values ​​and one maximum value (bistable and unstable). By tracking the trajectory of the current state point in the control variable space in real time, it is possible to intuitively see how the system moves from the safe zone to the danger zone.

[0060] In one implementation, the geometric distance from the current state point to the boundary of the bifurcation set of the nonlinear cusp catastrophe model is calculated based on the coordinates of the current state point in the control variable space. The instability prediction index is then calculated by combining the geometric distance and the energy entropy value, including the following steps: Based on the mathematical definition of the nonlinear cusp catastrophe model, a bifurcation set discrimination equation consisting of splitting factor control variables and normal factor control variables is constructed, and bifurcation set curves are plotted in the control variable space. An iterative search algorithm is used to calculate the shortest Euclidean distance from the coordinates of the current state point in the control variable space to the bifurcation set curve. After normalizing the energy entropy value, a weighting coefficient is assigned to the energy entropy value. The instability prediction index is calculated by combining the shortest Euclidean distance with the weighted energy entropy value and using the instability prediction index calculation formula.

[0061] In this implementation, based on the fundamental mathematical principles of cusp catastrophe theory, the system's potential energy function V(x) at the critical state has not only zero first derivative (equilibrium point) but also zero second derivative (inflection point), meaning the system has lost its restorative force to maintain its current state. By differentiating the standard potential energy function and solving the system of equations, eliminating the state variable x, we derive the algebraic relationship that the control variables u (split factor) and v (normality factor) must satisfy, i.e., the bifurcation set discrimination equation: The equation plots a characteristic curve on a two-dimensional control variable space plane with u as the horizontal axis and v as the vertical axis. This curve exhibits a sharp "V" shape or cusp shape and is called a bifurcation set curve. In the region inside the curve's opening, the system's potential function has two local minima, exhibiting a bistable state or an unstable jump region; while in the region outside the curve, the system has only one global minimum, exhibiting a stable monostable state.

[0062] After determining the bifurcation set curve as the stability boundary of the system, the current state point calculated in real time in the control variable space is... We need to find the shortest Euclidean distance from it to the bifurcation set curve. Since the bifurcation set curve is nonlinear, directly calculating the distance from a point to the curve analytically is difficult; therefore, a numerical iterative search algorithm is employed. The algorithm discretizes the bifurcation set curve, selecting a series of reference points, or constructs the objective function using gradient descent. Under constraints The algorithm searches for the minimum value. After rapid iterative convergence, the calculated shortest distance directly represents the stability margin of the system. The shorter the shortest distance, the closer the current state point is to the mutation boundary, the weaker the system's ability to resist disturbances and maintain stability, and the probability of a fault mutation increases exponentially.

[0063] Next, we need to introduce the energy entropy value of the friction texture. A comprehensive evaluation was completed. Since the entropy value may not be on the same order of magnitude as the geometric distance, it was first normalized. The historical maximum entropy value was then set. and minimum entropy value Using the linear transformation formula The entropy value is mapped to the [0,1] interval, making it a dimensionless scalar index. Then, based on practical engineering experience, a weighting coefficient is assigned to the normalized energy entropy value; this coefficient reflects the importance of microscopic wear characteristics in overall fault diagnosis. Finally, the features from the two dimensions are mathematically synthesized. The formula for calculating the instability prediction index is constructed as follows: In the formula, the first term uses the reciprocal of the geometric distance to characterize the risk of approaching the boundary; the smaller the distance, the larger the reciprocal. The first term is a very small constant to prevent the denominator from being zero; the second term is a weighted energy entropy value, with a larger entropy value indicating greater risk. The instability prediction index calculated by this formula integrates the vulnerability of the system's macroscopic dynamics (determined by geometric distance) and the degradation degree of the microscopic mechanical interface (determined by energy entropy). This index monotonically increases as the system state deteriorates. The index is low when the system is running smoothly, but rises sharply when the system approaches abrupt change boundaries or friction and wear intensify. This quantitative index provides the final decision-making basis for the control system, enabling fault prediction to move beyond fuzzy qualitative judgments and instead rely on a rigorous data-driven model, achieving high-precision real-time quantitative rating of the health status of the heave compensation device.

[0064] In one embodiment, when the instability prediction index meets a preset risk threshold condition, a variable damping control command is generated to adjust the hydraulic throttling element of the heave compensation device, and the heave compensation hydraulic cylinder of the heave compensation device is gradually locked to complete the fault self-protection, including the following steps: If the instability prediction index exceeds the preset risk threshold, the target damping coefficient is calculated using a fuzzy logic controller based on the extent to which the instability prediction index exceeds the risk threshold. The target damping coefficient is converted into an opening control signal for the hydraulic throttling element in the heave compensation device, and the proportional throttling valve of the heave compensation device is driven to reduce the flow area. During the adjustment process, the movement deceleration of the hydraulic cylinder of the heave compensation device is monitored in real time. If the movement deceleration exceeds the safety limit, the current throttling opening of the proportional throttle valve is maintained until the hydraulic cylinder completely stops moving and locks.

[0065] In this embodiment, when the real-time calculated instability prediction index exceeds the system's preset safety risk threshold, it indicates that the device has entered a potential fault zone, requiring a fuzzy logic control strategy for flexible intervention. First, the magnitude of the index exceeding the limit and the rate of change of that magnitude are calculated. These two variables are used as inputs to the fuzzy controller, mapped to linguistic variables such as "slight exceedance," "moderate exceedance," and "severe exceedance" through a predefined fuzzification interface. The fuzzy controller internally stores a set of inference rules based on expert experience, performs logical operations using Mamdani or Sugeno inference mechanisms, and outputs a fuzzy damping adjustment amount. Finally, defuzzification processing transforms the inference result into a precise numerical output, namely the target damping coefficient required under the current operating conditions. Based on the flow characteristic curve of the proportional throttle valve in the heave compensation device, an inverse mapping function between the damping coefficient and the valve opening is established. This function is typically non-linear, describing the specific position or flow area the valve spool needs to move to in order to generate a particular hydraulic resistance. The controller calculates the specific opening control signal based on this and sends it to the proportional throttle valve's drive circuit. Driven by current, the electromagnet generates thrust, overcoming the spring force to move the valve spool, thereby reducing the flow area of ​​the hydraulic oil through the valve port. As the flow area shrinks, the flow resistance in the hydraulic circuit increases, which in turn translates into a reverse damping force acting on the hydraulic cylinder piston rod.

[0066] While damping adjustment is being performed, the instantaneous deceleration of the hydraulic cylinder is monitored in real time using an accelerometer mounted on the cylinder or by second-order differentiation of displacement data. A safety limit based on the mechanical structural strength and personnel tolerance limits is set. During the process of controlling the proportional throttle valve to reduce its opening, the controller continuously compares the instantaneous deceleration with the safety limit. Once the current deceleration is detected to exceed the safety limit, it indicates that the damping force is increasing too rapidly, potentially causing water hammer in the pipeline or structural overload. The control logic immediately triggers a holding command, pausing the valve's further closing action, locking the current throttle opening, and waiting for the deceleration to fall back to a safe range before resuming adjustment. This dynamically constrained control strategy ensures that the entire shutdown process is smooth and controlled. As the speed gradually decreases, when the hydraulic cylinder speed is detected to eventually drop to zero or near zero, the system determines that a soft shutdown is complete, and then triggers a mechanical locking device or closes the main shut-off valve, safely locking the hydraulic cylinder in its current position, completing a shock-free fault self-protection process.

[0067] The present invention also discloses a fault prediction and self-protection system for diving bell heave compensation, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the fault prediction and self-protection method for diving bell heave compensation as described in any of the above.

[0068] The processor can be a central processing unit (CPU). Of course, depending on the actual use, it can also be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), off-the-shelf programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor can be a microprocessor or any conventional processor, etc., and this application does not limit it.

[0069] The memory can be an internal storage unit of a computer device, such as a hard disk or RAM, or an external storage device, such as a plug-in hard disk, smart memory card (SMC), secure digital card (SD), or flash memory card (FC) provided on the computer device. Furthermore, the memory can be a combination of internal storage units and external storage devices of a computer device. The memory is used to store computer programs and other programs and data required by the computer device. The memory can also be used to temporarily store data that has been output or will be output. This application does not limit this.

[0070] The present invention also discloses a computer-readable storage medium storing instructions that, when executed by a processor, configure the processor to perform the fault prediction and self-protection method for diving bell heave compensation described in any of the above embodiments.

[0071] The computer program can be stored in a machine-readable medium. The computer program includes computer program code, which can be in the form of source code, object code, executable file, or certain middleware. The machine-readable medium includes any entity or device capable of carrying computer program code, recording media, USB flash drive, portable hard drive, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signals, telecommunication signals, and software distribution media, etc. It should be noted that the machine-readable medium includes, but is not limited to, the above-mentioned components.

[0072] The fault prediction and self-protection method for diving bell heave compensation in the above embodiments is stored in the computer-readable storage medium and loaded and executed on the processor to facilitate the storage and application of the above method.

[0073] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of protection of this application is limited to these examples; within the framework of this application, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of one or more embodiments of this application as described above, which are not provided in detail for the sake of brevity.

[0074] One or more embodiments in this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of this application. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of one or more embodiments in this application should be included within the protection scope of this application.

Claims

1. A method for failure prediction and self-protection of a diving bell heave compensation, characterized in that, Includes the following steps: Acquire real-time operating data of the heave compensation device, extract hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle from the real-time operating data, map them to a two-dimensional coordinate system to construct a closed hysteresis loop, and convert the closed hysteresis loop into a binary hysteresis feature field image. Geometric moment analysis is performed on the hysteresis feature field image, and topological feature descriptors for quantizing the discreteness and skewness characteristics of the hysteresis loop are calculated. Homomorphic deconvolution is used to process the accumulator pressure data and transform it to the cepstrum domain. Low-frequency components of the accumulator pressure data are filtered out by inverted high-pass filtering. At the same time, the friction texture signal characterizing the mechanical contact state is separated and the energy entropy value of the friction texture signal is calculated. The discreteness characteristics are mapped to splitting factor control variables, and the skewness characteristics are mapped to normal factor control variables. A control variable space consisting of splitting factor control variables and normal factor control variables is constructed, and the current state point of the heave compensation device in the control variable space is determined. The geometric distance from the current state point to the boundary of the bifurcation set of the nonlinear cusp catastrophe model is calculated based on the coordinates of the current state point in the control variable space. The instability prediction index is then calculated by combining the geometric distance and the energy entropy value. When the instability prediction index meets the preset risk threshold condition, a variable damping control command is generated to adjust the hydraulic throttling element of the heave compensation device, and the heave compensation hydraulic cylinder of the heave compensation device is gradually locked to complete the fault self-protection.

2. The method for failure prediction and self-protection of a bell heave compensation according to claim 1, characterized in that, The step of extracting hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle from real-time operational data, mapping them to a two-dimensional coordinate system to construct a closed hysteresis loop, and converting the closed hysteresis loop into a binary hysteresis feature field image includes the following steps: Set a two-dimensional grid resolution that is adapted to the displacement range of the hydraulic cylinder and the pressure variation range of the accumulator, and initialize the two-dimensional grid. The two-dimensional grid resolution is used to define the pixel density of the hysteresis feature field image. The hydraulic cylinder displacement data and accumulator pressure data within a single heave cycle in the real-time operation data are processed by cubic spline interpolation to generate a smooth and continuous closed hysteresis loop trajectory. The coordinates of the closed hysteresis loop trajectory are normalized to a unit interval, and the normalized coordinates are mapped to a two-dimensional grid. The grid cells passed by the closed hysteresis loop trajectory are marked as foreground pixels, and the grid cells not passed by the closed hysteresis loop trajectory are marked as background pixels, thus generating an initial two-dimensional grid image. The initial two-dimensional mesh image is preprocessed using morphological closing operations to obtain a two-dimensional mesh image. The two-dimensional mesh image is then transformed into a binary matrix containing only 0 and 1 values ​​to generate a hysteresis feature field image for feature extraction.

3. The method for the fault prediction and self-protection of the diving bell heave compensation according to claim 2, characterized in that, The step of performing geometric moment analysis on the hysteresis feature field image and calculating the topological feature descriptors for quantizing the discreteness and skewness features of the hysteresis loop includes the following steps: The zeroth moment of the hysteresis feature field image is calculated to obtain the closed area of ​​the hysteresis loop, and the geometric centroid coordinates of the hysteresis loop in the two-dimensional coordinate system are determined by calculating the first moment of the hysteresis feature field image. The second and third central moments of the hysteresis feature field image are calculated based on the geometric centroid coordinates. The second and third central moments are then normalized using the zero-order moment to obtain the normalized central moments. The first invariant characteristic formula is constructed using the second-order central moment in the normalized central moments, and the sum of the moments of inertia in the orthogonal directions is defined as the discreteness characteristic that characterizes the discreteness of the figure. The second invariant feature formula is constructed using the third central moment in the normalized central moments. The sum of the projection skewness in the orthogonal directions is defined as the skewness feature that characterizes the degree of skewness of the figure. The calculated discreteness feature and the skewness feature are combined to form a topological feature descriptor.

4. The bell heave compensation fault prediction and self-protection method according to claim 1, characterized in that, The process of processing the accumulator pressure data using homomorphic deconvolution and transforming it to the cepstrum domain, filtering out low-frequency components of the accumulator pressure data using an inverted high-pass filter, and simultaneously separating the friction texture signal characterizing the mechanical contact state and calculating the energy entropy value of the friction texture signal includes the following steps: A Hamming window is applied to the collected accumulator pressure data, and a fast Fourier transform is performed on the windowed accumulator pressure data to obtain the pressure spectrum. Taking the natural logarithm of the amplitude of the pressure spectrum transforms the convolutional relationship in the time domain of the pressure spectrum into an additive relationship in the logarithmic amplitude spectrum in the frequency domain. Perform an inverse fast Fourier transform on the logarithmic amplitude spectrum in the frequency domain to map the accumulator pressure data to the inverted frequency domain to obtain the cepstrum sequence of the accumulator pressure data; The inverted frequency cutoff threshold is set according to the wave compensation period of the heave compensation device, and an inverted high-pass filter is constructed to filter out the low inverted frequency components representing wave motion and fluid pressure in the cepstral sequence. An inverse homomorphic transform is performed on the high inverse frequency residual component after processing by the inverse high-pass filter, and the high inverse frequency residual component is restored to the time domain to obtain the friction texture signal; The probability density function of the friction texture signal is calculated, and the entropy value of the probability density function is calculated based on the Shannon entropy formula to obtain the energy entropy value that quantifies the complexity of the mechanical contact surface.

5. The bell heave compensation fault prediction and self-protection method according to claim 4, characterized in that, The step of setting the inverted frequency cutoff threshold according to the wave compensation period of the heave compensation device and constructing an inverted high-pass filter to filter out the low inverted frequency components representing wave motion and fluid pressure in the cepstral sequence includes the following steps: Calculate the amplitude envelope of the cepstrum sequence and normalize the amplitude envelope to construct a histogram of the cepstrum energy probability distribution; The information complexity of the distribution histogram is quantified by Raney entropy, and the sum of the Raney entropy of the background class cepstral distribution and the foreground class cepstral distribution is defined as the objective function of the candidate cepstral cutoff threshold. Traverse all candidate cepstral frequency cutoff thresholds within the domain of the cepstral sequence, solve the global maximum value of the objective function using a genetic algorithm, and determine the candidate cepstral frequency cutoff threshold corresponding to the global maximum value as the optimal cepstral frequency cutoff threshold; Based on the optimal cepstral frequency cutoff threshold, a cosine tapered window function with smooth transition characteristics is constructed. The cosine tapered window function is multiplied with the cepstral sequence to suppress wave fluid interference at low cepstral frequencies while preserving the friction characteristics of high cepstral frequencies.

6. The fault prediction and self-protection method for diving bell heave compensation according to claim 1, characterized in that, The steps of mapping discreteness features to splitting factor control variables, mapping skewness features to normality factor control variables, constructing a control variable space composed of splitting factor control variables and normality factor control variables, and determining the current state point of the heave compensation device in the control variable space include the following: A nonlinear cusp mutation model of the heave compensation device is constructed and a standard cusp mutation potential function is established. The standard cusp mutation potential function includes a state variable, a quartic term, a quadratic term weighted by a splitting factor, and a linear term weighted by a normality factor. Historical operating data of the heave compensation device under healthy conditions were collected, and a baseline distribution model was established by calculating the dispersion and skewness characteristics of the historical operating data. Calculate the first standardized Z-score of the discreteness feature extracted at the current time relative to the baseline distribution model, and the second standardized Z-score of the skewness feature extracted at the current time relative to the baseline distribution model; The first standardized Z-score is compressed and mapped using a nonlinear Sigmoid function to generate a splitting factor control variable corresponding to the quadratic term of the standard cusp mutation potential function. The hyperbolic tangent function is used to compress and map the second standardized Z-score, generating a normal factor control variable corresponding to the first term of the standard cusp mutation potential energy function; A control variable space consisting of splitting factor control variables and normal factor control variables is constructed, and the current state point of the heave compensation device in the control variable space is determined based on the standard cusp mutation potential energy function.

7. The fault prediction and self-protection method for diving bell heave compensation according to claim 6, characterized in that, The process of calculating the geometric distance from the current state point to the boundary of the bifurcation set of the nonlinear cusp catastrophe model based on the coordinates of the current state point in the control variable space, and then combining the geometric distance and the energy entropy value to calculate the instability prediction index includes the following steps: Based on the mathematical definition of the nonlinear cusp catastrophe model, a bifurcation set discrimination equation consisting of splitting factor control variables and normal factor control variables is constructed, and bifurcation set curves are plotted in the control variable space. An iterative search algorithm is used to calculate the shortest Euclidean distance from the coordinates of the current state point in the control variable space to the bifurcation set curve. After normalizing the energy entropy value, a weighting coefficient is assigned to the energy entropy value. The instability prediction index is calculated by combining the shortest Euclidean distance with the weighted energy entropy value and using the instability prediction index calculation formula.

8. The fault prediction and self-protection method for diving bell heave compensation according to claim 7, characterized in that, When the instability prediction index meets the preset risk threshold condition, a variable damping control command is generated to adjust the hydraulic throttling element of the heave compensation device, and the heave compensation hydraulic cylinder of the heave compensation device is gradually locked to complete the fault self-protection, including the following steps: If the instability prediction index exceeds the preset risk threshold, the target damping coefficient is calculated using a fuzzy logic controller based on the extent to which the instability prediction index exceeds the risk threshold. The target damping coefficient is converted into an opening control signal for the hydraulic throttling element in the heave compensation device, and the proportional throttling valve of the heave compensation device is driven to reduce the flow area. During the adjustment process, the movement deceleration of the hydraulic cylinder of the heave compensation device is monitored in real time. If the movement deceleration exceeds the safety limit, the current throttling opening of the proportional throttle valve is maintained until the hydraulic cylinder completely stops moving and locks.

9. A fault prediction and self-protection system for diving bell heave compensation, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the fault prediction and self-protection method for diving bell heave compensation as described in any one of claims 1 to 8.

10. A computer-readable storage medium storing instructions thereon, characterized in that, When executed by the processor, the instruction causes the processor to be configured to perform the fault prediction and self-protection method for diving bell heave compensation according to any one of claims 1 to 8.

Images