Modal Parameters Identification of A Real Offshore Platform From the Response Excited by Natural Ice Loading

This paper investigates the possibility of utilizing response from natural ice loading for modal parameter identification of real offshore platforms. The test platform is the JZ20-2MUQ jacket platform located in the Liaodong Bay, China. A field experiment is carried out in winter season, as the platform is excited by floating ices. The feasibility is demonstrated by the acceleration response of two different segments. By the SSI-data method, the modal frequencies and damping ratios of four structural modes can be successfully identified from both segments. The estimated information from both segments is almost identical, which demonstrates that the modal identification is trustworthy. Furthermore, by taking the Jacket platform as a benchmark, the numerical performance of five popular time-domain EMA methods is systematically compared from different viewpoints. The comparisons are categorized as: (1) stochastic methods versus deterministic methods; (2) high-order methods versus low-order methods; (3) data-driven versus covariance-driven stochastic subspace identification methods.


Introduction
Accurate knowledge of the modal parameters is crucial to many practical engineering issues, for example vibration control, finite element model updating and structural health monitoring. To extract modal properties of a structure from measured vibration data has been termed the experimental modal analysis (EMA) (Maia and Silva, 1997;Ewins, 2000). EMA, initiated from early 1970's, has obtained significant progress and been successfully applied to mechanical, automotive, civil and aerospace engineering (Zhang, 2013). Classical EMA techniques are input-output methods, which have been developed based on the measurement of vibration response, together with the deterministic knowledge of the artificial input excitation (Cunha and Caetano, 2006). In contrast, without the deterministic knowledge of the input, modern output-only modal identification methods, developed from 1990's, would utilize response data only. They are also named as operational modal analysis (OMA) methods. OMA means that the modal testing and analysis for a structure is in its operational condition subjected to natural excitation with output-only measurements (Peeters and de Reock, 2001;Yang et al., 2015). Although the deterministic input knowledge is not required for output-only methods, they still need to have a mathematical assumption about the stochastic properties of the excitation. All current output-only modal identification methods assume that the input is a zero-mean Gaussian white noise random process (Peeters and de Roeck, 1999). It should be pointed out that, in practical application, these mathematical assumptions are usually not rigorously obeyed. In contrast, it has been widely accepted that stochastic modal analysis methods may have the ability to identify modal parameters from a deterministic case, and vice versa (Peeters and Ventura, 2003). In the past decades, a number of EMA/OMA methods have been developed and some of them have been applied to offshore platforms (Brincker et al., 1996;Mangal et al., 2001;Yang et al., 2006;Hu et al., 2011;Li et al., 2012;Liu et al., 2015).
Regarding the excitation for modal testing, an offshore platform could be impacted by a maneuvered object such as a small boat or be pulled and then released after a static displacement is reached (Mangal et al., 2001;Li et al., 2012). As the size of offshore platforms becomes huger, imposing artificial excitation is usually time-consuming, expensive, sophisticated and even impossible. In this situation, it is suggested to employ the response from natural excitation for modal analysis. Modal parameters of some platforms have been extracted by utilizing the responses from wind loading or wave loading (Brincker et al., 1996;Yang et al., 2006;Liu et al., 2015). However, to employ the response from natural ice loading for modal parameter identification is rarely reported, though exploiting oil and gas in cold regions is attracting more attention from the offshore engineering industry. For example, in Canadian and Alaskan Beaufort Sea, the Barents and Kara Seas, the Liaodong Bay, and offshore Sakhalin, ices could be widespread (Gerwick, 2007). Venturella et al. have employed the modal analysis technique for the ice-structure interaction problem, but their work is a forward problem of analytical modal analysis, in which the response of each mode is superposed to find the complete modal response of the entire length of a pier subject to incremental ice loading (Venturella, 2008;Venturella et al., 2008b).
Regarding the methods for modal analysis, modal parameters of offshore platforms have been probed by many techniques, ranging from the basic Fourier-based peak-picking method to the sophisticated stochastic subspace identification method (Xin et al., 2012). However, a comparative study on the performance of different EMA/OMA methods specific for offshore platforms is still lacking. Meanwhile, this kind of comparative studies have been carried out for other engineering structures. For instance, Abdelghani et al. (1998) presented a detailed comparison study of three subspace identification methods on the NASA Mini-Mast structure. Florakis et al. (1998) addressed a critical comparison and assessment of a deterministic ARX, a stochastic ARMAX, and the deterministic ERA structural dynamics identification methods based on vibration data obtained from a scale airplane structure. Petsounis and Fassois (2001) conducted a critical comparison of four stochastic and three deterministic methods for the parametric time-domain identification of a railway structure from random excitation and noise-corrupted response signals. Peeters and Ventura (2003) evaluated the performance of different modal analysis methods for analyzing the dynamic characteristics of a three-span reinforced concrete bridge in Switzerland from forced, free and ambient vibration data.
There are two objectives in this paper. The primary objective is to investigate the applicability and the reliability of utilizing responses from natural ice loading for modal parameter identification of offshore platforms. Another objective is to compare the numerical performance of popular time-domain modal analysis methods on a real offshore platform. Specifically examined are the poly-reference complex exponential (PRCE) method, the eigensystem realization algorithm (ERA) method, the instrumental variable (IV) method, the covariance-driven stochastic subspace identification (SSI-cov) method and the data-driven stochastic subspace identification (SSI-data) method. The PRCE and ERA are deterministic methods while the other three are stochastic (Peeters and de Reock, 2001). Hu et al. (2014) revealed that, as various timedomain methods are modeled differently with distinct solution algorithms, implementing these methods has different performance because of the numerical conditioning and stability. However, the numerical observation was based on a case study for a simulated system of 6 degree-of-freedom. It is of interest to further compare the performance of different methods on real structures.
The platform to be tested is the JZ20-2MUQ jacket platform located in the Liaodong Bay, China. A field experiment is carried out in winter season, when the platform is excited by floating ices. An accelerometer is installed on a leg and the acceleration response is recorded. In this paper, the modal parameters of the platform will be identified from two different acceleration segments by adopting the SSIdata method. For both acceleration segments, the response is measured immediately after the platform is impacted by a floating ice. The purpose of utilizing two different segments is to validate the identification results from each other. Secondly, the numerical performance of the above-mentioned methods will be compared from different viewpoints, including deterministic methods versus stochastic methods, high-order methods versus low-order methods, and covariance-driven stochastic subspace method versus data-driven stochastic subspace method. To distinguish the structural modal information, stability diagrams are constructed by identifying parametric models of increasing order.

Models for modal identification methods
Both deterministic methods and stochastic methods can be divided as parametric methods essentially developed in the time domain and nonparametric methods in the frequency domain. The frequency-domain methods always suffer the leakage errors in the estimates and the frequency resolution issue of spectral estimates. Employing time-domain parametric methods could provide better accuracy for modal parameter estimation, especially when a large frequency range or a large number of modes exist in the data. For each time domain EMA method, the first step is to adopt an appropriate mathematical model to represent the dynamic system. This section reviews four mathematical models, of which two are deterministic and the other two stochastic. The models are based on either a high-order matrix polynomial difference equation or a first-order state-space equation.

Deterministic models
In deterministic models, the inputs are known without any stochastic assumptions. For the PRCE method, the corresponding deterministic model is : Mathematically, it is a high-order matrix finite difference equation, where , , and are the model order, the number of time steps, the j-th matrix coefficient, and the impulse response function or free relaxation response matrix, respectively. The deterministic model associated with the ERA method is of first-order (Juang, 1994): It is a state-space model. In this model, and represent the state vector and the state matrix; and are the input vector and the input influence matrix; and and are the output vector and output influence matrix.

Stochastic models
Output-only modal analysis is usually employed when the input is unknown but can be assumed as a white noise random process. Similar to the deterministic models, there are also a high-order stochastic model and a first-order stochastic state-space model. The high-order model is formulated as (Ljung, 1999): It is commonly known as the ARMA model, with the Auto-Regressive (AR) part at the left-hand side and the Moving-Average (MA) part at the right-hand side of the equation. The matrices and are the AR matrix parameters and the MA matrix parameters, respectively. and represent the output vector and a white noise vector sequence, respectively.
As the counterpart of Eq.
(2), when the system is excited by a white noise random process, the stochastic statespace model can be written as (Peeters and de Reock, 2001): where and are the process noise and the measurement noise, respectively. Since they are unmeasurable, both of them are assumed to be zero-mean white process. The covariance matrices between them can be expressed as: where is the expected value operator and is the Kro-necker delta.

Algorithms for modal identification methods
This section summarizes the algorithms of the methods investigated and compared in this paper, which include the poly-reference complex exponential (PRCE) method, the eigensystem realization algorithm (ERA) method, the instrumental variable (IV) method, the covariance-driven stochastic subspace identification (SSI-cov) method and the data-driven stochastic subspace identification (SSI-data) method.

Polyreference Complex Exponential (PRCE) method
PRCE method is based on the high-order matrix polynomial model (see Eq. (1)) and there are two stages in the algorithm. The first stage is to calculate the 2N unknown poles from the model. For simplification, the identity matrix I is chosen for the matrix coefficient of the highest order, b n . The matrix form for the equations in Eq. (1) can be expressed as (Maia and Silva, 1997;: In real application, the system is usually over-determined, which means , where is the order of the matrix polynomial. The coefficient matrices can be solved by the least squares method in the over-determined cases. Once the coefficient matrices are known, the system poles could be obtained by computing the roots of the characteristic polynomial equation where , being the system eigenvalues and Δt, the sampling interval. The companion matrix approach is usually applied to obtain these roots (Maia and Silva, 1997;Allemang and Brown, 1998).

Eigensystem Realization Algorithm (ERA) method
When a unit impulse excitation is taken for each input element, the results from Eq. (2) can be assembled into a pulse-response matrix as follows (Juang, 1994): The constant matrices Y k in the sequence are the transpose of h k and termed as Markov parameters. The notation Y k is employed following the origin of the ERA method (Juang, 1994).
For the ERA method, a system realization is the computation of a triplet [A, B, C] from the Markov parameters shown in Eq. (11), for which the discrete-time model, Eq.
(2), is satisfied. Any system has an infinite number of realizations which will predict the identical response for any particular input, and a minimum realization means a model with the smallest state-space dimensions among all realizable systems that have the same input-output relations.
The first step in the system realization is to form the generalized block Hankel matrix, with the Markov parameters from Eq. (11): ξ η where and are the number of block rows and columns in the block Hankel matrix, respectively. Differing from other system realization methods which use the generalized Hankel matrix given in Eq. (12), the ERA algorithm begins by forming an ERA block data matrix which is obtained by deleting some rows and columns of the generalized Hankel matrix of Eq. (12), but maintaining the first block matrix, Y k , intact. Furthermore, the standard ordering of entries in the generalized Hankel matrix does not need to be maintained (Juang, 1994). However, for the simplicity of presentation in this paper, the block Hankel matrix given in Eq. (12) is used as the ERA block data matrix.
After substituting the Markov parameters from Eq. (11) into Eq. (12), H(k-1) can be decomposed into three matrices: where the block matrix is the observability matrix: Q η and the block matrix is the controllability matrix: In particular, substituting k = 1 and k = 2 into Eq. (13) yields and (17) The ERA process starts with the factorization of H(0), which is obtained by replacing k = 1 in Eq. (12), using singular value decomposition where the columns of matrices U and V are orthonormal and U 1 and V 1 are the matrices formed by the first n columns of U and V, respectively; with monotonically non-increasing , i.e. . Comparison of Eqs. (16) and (18) suggests that is related to U 1 and is related to . Indeed, one possible choice is and . This choice appears to make both and balanced. With substituting these choices for and into Eq. (17), one can obtain that A Q η P ξ where the quantity means the estimated A to distinguish from the true quantity. From Eqs. (14) and (15), it is clear that the first N i columns of form the input matrix B whereas the first N o rows of form the output matrix C.

Instrumental Variable (IV) method
For the same structural system, the response from the impulse loading and the output covariances from a white noise excitation have the similar mathematical expressions de Roeck, 1999, 2001). This principle has helped to extend the classical experimental modal analysis methods which were based on impulse response to utilize the output covariances.
The IV method is also based on the ARMA model. It assumes that the residuals e k are uncorrelated with past data and then extracts the maximum data from the response. It estimates only the parameters for the AR part and the equations of IV method can be expressed in terms of the output covariances as (Ljung, 1999): One can clearly see that this equation and the mathematical model of the PRCE method are similar, and the only difference is that PRCE utilizes impulse response and IV employs output covariances ( ). Hence, they can apply the same algorithm for calculating the system eigenvalues.

Covariance-driven Stochastic Subspace Identification
(SSI-cov) method Referring to Eq. (4), one assumes the response stochastic process to be stationary with and , where the state covariance matrix is independent of the time . This implies that is a stable matrix. Since and are zero mean white noise vector sequences, independent of , one has and . Then one finds the Lyapunov equation for the state covariance matrix: . Defining the covariance matrix between and as , one can obtain: , and The classical algorithm of the SSI-cov method is to arrange the output covariance matrices R i in a Toeplitz matrix form: is the extended observability matrix and is the reversed extended stochastic controllability matrix. They can be written as: There are several ways to get a minimum realization for the system matrix A. One way is to use only the extended observability matrix (see Eq. (24)). Let the singular value decomposition of to be expressed as: In practice, the order of the system, n, is set equal to the number of singular values in Eq. (25) above a preset threshold. One possible realization of is . The matrix A can now be determined from the extended observability matrix by making use of the shift structure of the matrix . Denote as the extended observability matrix without the first l rows, and the extended observability matrix without the last l rows. Observing Eq. (24), one can easily find , which suggests that an estimate for A can be computed as: where the superscript † denotes the Moore-Penrose pseudoinverse of a matrix. The estimated C could be obtained from the first block of . Because flipping left to right a Toeplitz matrix (see Eq. (23)) yields a Hankel matrix, a realization of A can also be estimated from the information of covariance matrices by adopting the same steps of the Eigensystem Realization Algorithm (ERA) (Peeters and de Reock, 2001).

Data-driven Stochastic Subspace Identification (SSIdata) method
The data-driven stochastic subspace identification (SSIdata) method avoids the computation of covariances between the outputs. It is replaced by projecting the row space of future outputs into the row space of past outputs. A detailed derivation of its algorithm is given in van Overs-chee and de Moor (1996).
The output measurements , are gathered in a block Hankel matrix H 0|2i−1 with 2i block rows and j columns, where j = s − 2i + 2: The matrices Y p and Y f are defined by splitting H 0|2i−1 into two parts of i block rows, in which the subscripts p and f stand for past and future. The projection of the row space of the future outputs Y f into the row space of the past outputs Y p , denoted by , can be computed by: ℘ i However, a more accurate way to obtain numerically is via the QR-factorization of H 0|2i−1 , instead of computing the above formula directly.
In the theoretical derivations of the stochastic subspace identification, one assumes that the number of measurements goes to infinity. Assuming j → ∞, the main theorem of stochastic subspace identification states that the projection is equal to the product of the observability matrix and the Kalman filter state sequence (van Overschee and De Moor, 1996): where the forward Kalman filter state sequence is .
This theorem can algebraically be summarized as follows: (1) rank of = n, (2) row space of = row space of , and (3) column space of = column space . This summary is the essence of why the projection-based algorithm has been called a subspace algorithm, as it retrieves system related matrices as subspaces of projected data matrices. Several approaches can be utilized to estimate A and C. One of them is to use the extended observability matrix alone (see Eq. (26)). Another approach uses the Kalman filter state sequence together with the extended observability matrix (van Overschee and de Moor, 1996).

Identification of true modes for structures
The ultimate goal of the experimental modal analysis is to correctly identify "true modes" and accurately estimate the corresponding modal parameters. As described above, various methods identify the modal parameters through different models and distinct solution algorithms. The application of the methods would have different numerical conditioning and numerical stability, which are relevant to the perturbation behavior of a mathematical problem (model) itself and the perturbation behavior of an algorithm used to solve that problem on a computer, respectively (Trefethen and Bau, 1997). As a result, these modal identification methods must have different numerical performances on extracting the true structural modes.
For identifying true modes, one often relies on using stability diagrams. The stability diagrams involve tracking the estimates of frequency, damping, and possibly other modal properties as a function of model order. The poles corresponding to a certain model order are compared with the poles of a one-order-lower model. If the modal frequency, the damping ratio and the related mode shape differences are within preset limits, the pole is labeled as a stable one. As the model order is increased, more and more modal frequencies are estimated but, hopefully, the estimates of the physical modal parameters will stabilize. For modes that are very active in the measured data, the modal parameters will stabilize at a very low model order. For modes that were poorly excited in the measured data, the modal parameters may not stabilize until a very high model order is chosen. Nevertheless, the traditional wisdom has suggested that the spurious "computational" poles will not stabilize at all during this process.

Field test and modal parameter identification
The test structure is the JZ20-2MUQ platform, a fourlegged jacket-type offshore platform located in the northern region of the Liaodong Bay, China (see Fig. 1), which was installed in Aug., 1992. Its plane dimensions are approximately 25 m×24 m at the deck elevation, which is 15.5 m above the still water level. With reference to the still water level, the elevations of the structure's three stories are +5.85 m, −3.5 m, and −15.5 m, respectively. The four sloping legs, which have an inclination of 10:1, are cylindrical piles with a uniform outer diameter of D=1.52 m and a wall thickness of t=0.025 m. Because of the high latitude, it may endure the action of ice loading for 2 or 3 months every year. A field test was carried out. In the experiment, a tri-axial accelerometer was installed on a leg at the vertical position of +10 m, recording the vibration signals in horizontal and vertical directions, respectively. The position of the accelerometer is shown as a red circle in Fig. 1.
Two segments of response from a horizontal direction (see Fig. 2), one of 30-second and the other one of 40second with the sampling frequency of 500 Hz, are utilized for modal parameter identification in this paper. Each segment is chosen after the platform is impacted by a floating ice. During that process, the excitation to the platform can be assumed as a combination of the impact force and a white noise loading. Meanwhile, in the short time of excitation, the loading can be assumed quasi-linear. The excitation basically satisfies the requirement of the numerical models and algorithms of modal parameter identification methods. The purpose of employing two segments is to val-idate the identification results from each other. This crossvalidation has much significance for the test platform. As the test platform was built more than twenty years ago, its modal parameters might have changed from its original state due to various reasons, such as marine growth, structure deterioration, etc. Consequently, the "analytic" modal parameters are unknown.
In this section, one of the most advanced methods, the data-driven stochastic subspace identification (SSI-data) method is applied to identify the modal parameters. Part of the identification results has been reported by the authors previously (Yang et al., 2016). The excellent numerical performance of SSI-data method originates from its first-order stochastic model and the QR algorithm for computing the structural eigenvalues (Boonyapinyo and Janesupasaeree, 2010;Hu et al., 2014). Although this method is mathematically derived based on the assumption that the structural system is under the excitation of Gaussian white noise, the theoretical assumption of white noise turns out to be not too strict in practical applications. It has been demonstrated that as long as the unknown input spectrum is quite flat, this method can accurately identify the modal parameters of a structure (Peeters and Ventura, 2003).  Traditionally, the first step to employ a time-domain modal analysis method is determining the model order. For SSI-data method, it could be achieved by analyzing the normalized singular value of the projection matrix . A conventional way to choose the model order is to find a significant gap of the normalized singular value. In this study, the Hankel matrix for the SSI-data method has been constructed with i=256. The selection of this size is not arbitrary, but a careful choice based on multiple trials with different sizes of the Hankel matrix. It was found that the selection of i=256 could identify more modes with better estimated values. When i=256, the curves of the normalized singular values of the projection matrix for the two segments are shown in Fig. 3. Observing this figure, the singular values approach to an asymptote after the model order gets 8 for both segments, which means no more than 4 modes can be identified from the signals. By taking model order as 8, the estimated modal frequencies from segment 1 are 0.91 Hz, 3.06 Hz, 10.81 Hz and 22.47 Hz. The corresponding damping ratios are 3.29%, 2.83%, 1.42% and 2.48%, respectively. For segment 2, the estimated modal frequencies are 0.89 Hz, 3.07 Hz, 21.08 Hz and 48.65 Hz. The corresponding damping ratios are 2.01%, 1.97%, 6.57% and 1.78%, respectively. It is obvious that three modal frequencies are in high agreement, but there are some discrepancies associated with other parameter values. It is clear that the identified mode with the frequency of 10.81 Hz is a fake mode from the noise in the responses. There are mainly two ways to eliminate the influence of the noise. One way is to remove the noise from the response before using it for modal parameter identification Bao and Shi, 2019;Bao et al., 2020). Another way is to use stability diagram.
Nowadays, one is less interested in a good model order determination, but rather in the modal parameters extracted from running models with different orders. Past experience with real data showed that it was better to overspecify the model order and to eliminate spurious numerical poles afterwards. This could be done by constructing stability diagrams. When producing a stability diagram, the poles cor- responding to a certain model order are compared with the poles of a one order-lower model. If the modal frequency and the damping ratio differences are within preset limits, the pole is labeled as a stable one. For the stability diagrams in this paper, the poles are labeled as stable when they are within the limitation of 1% difference for frequency and 5% for damping ratio among two consecutive model orders, namely and , where and denote an estimated frequency and damping ratio from model order n. For all stability diagrams throughout this paper, the poles with symbol '*' are for stable ones, and '•' for unstable ones. The background curve of the stabilization diagram is a Fourier spectrum, which is obtained from the DFT (discrete Fourier transform) of the sample signal.
With the model order varying from 2 to 40, the stability diagrams obtained from implementing the SSI-data method on both segments are shown in Fig. 4 and Fig. 5, respectively (Yang et al., 2016). A glimpse on these figures tells one that the stability diagrams from the SSI-data method are very clean. From the stability diagrams, four modes can be easily identified. Besides, the poles associated with the modes can get stabilized at low model order, which rein- Fig. 3. Singular value curves associated with the projection matrix .  forces the traditional wisdom that the SSI-data method is a powerful tool for modal analysis. For more information, the values of estimated frequencies and damping ratios related to the four identifiable modes at different model orders are listed in Table 1 and Table 2 (Yang et al., 2016). In these tables, "N/A" means the corresponding mode is not available and "N/S" means the associated poles are not stable. From the tables, one can see that the estimated values from two segments are very close to each other, especially the modal frequencies. By taking the average value from the stable poles, the modal frequencies from Segment 1 can be determined as 0.91 Hz, 3.06 Hz, 7.67 Hz and 22.57 Hz. The corresponding damping ratios can be determined as 3.60%, 2.99%, 2.90% and 2.30%, respectively. For Segment 2, the modal frequencies can be determined as 0.89 Hz, 3.07 Hz, 7.67 Hz and 22.53 Hz. The corresponding damping ratios can be determined as 3.21%, 2.75%, 2.71%, 2.82%, respectively. The estimated frequencies of four identifiable modes from each segment are almost identical, which can demonstrate that the modal parameter identification results are reliable. The discrepancies among the identified damping ratios are within the expectation, as it has been recognized for a long time that the identification results of damping ratios are of low accuracy (Ewins, 2000;Peeters and Ventura, 2003).
It must be pointed out that not all the data from the ac-celerometer are perfect for modal identification. Multiple segments of the recorded signal are tried and from the trial the above-mentioned two segments are analyzed detailedly in this paper. The effectiveness of the data depends on the ice condition and the attribute of the impact loading. Based on the authors' experience, it is desired to employ the measurement data after the leg of the platform is impacted by a floating ice with high momentum and the impact is finished in a short time.

Comparison on the numerical performance of EMA methods
In the previous section, the applicability of utilizing response from natural ice loading is verified by the SSI-data method and two different acceleration response segments. In this section, the numerical performance of the PRCE method, the ERA method, the IV method, the SSI-cov method and the SSI-data method will be compared from different viewpoints: (1) The deterministic methods are compared with their counterparts in stochastic modal analysis, which means the PRCE method is compared with the IV method and the ERA method is compared with the SSI-cov method.
(2) The methods based on high-order model are compared with the methods based on first-order model. Specifically, the PRCE method is compared with the ERA method for deterministic modal analysis, and the IV method is compared with the SSI-cov method for stochastic modal analys- is.
(3) The two stochastic subspace identification (SSI) methods are compared, i.e. the SSI-cov method is compared with the SSI-data method. All of these methods have been tested on both Segment 1 and Segment 2 and it is found that their performances on two segments are quite similar. To be concise, only the identification results from Segment 1 are presented in details in this section.

Deterministic methods versus stochastic methods
It is well-known for some time that there exist similar mathematical expressions for deterministic and stochastic methods. Although derived in a different way, the final equations of the IV method and the SSI-cov method correspond to the PRCE method and the ERA method, respectively. In modal analysis applications, the distinction is to utilize the response directly or the output covariance instead.
C i With the model order ranging from 2 to 40, the stability diagrams from the PRCE method, the IV method, the ERA method and the SSI-cov method are shown in Figs. 6−9. It should be mentioned that these stability diagrams are chosen after the parameter study of each method has been fulfilled. The purpose is to find the best estimation result from each method and then compare them fairly and objectively. While constructing the stability by the ERA method, different sizes of Hankel matrix in Eq. (12) are tried. Shown in Fig. 8 is the stability diagram with the number of rows in Hankel matrix as 100, which has a relatively better performance than others. In the computation of the stability diagrams by the IV and the SSI-cov methods, only the first 256 time steps of the computed covariance function are used for the estimation of the modal parameters. That is because the covariance values with a smaller time lag have better accuracy when the measured response is of finite length. In the SSI-cov method, when 256 covariance values are employed, the number of rows of Toeplitz matrix is 128 (i = 128).
Comparison between Figs. 6 and 7 suggests that IV method performs obviously better than the PRCE method.
From Fig. 7, the modal frequencies of 3.06 Hz and 22.57 Hz can be easily identified by IV method. In contrast, stable poles from PRCE method for the frequency of 3.06 Hz are not formed even when the model order gets 40 (see Fig. 6). Though the PRCE method can also roughly identify the frequency of 22.57 Hz, the corresponding poles do not become stable until the model order equals 30. Besides, the stability diagram from IV method (see Fig. 7) is a little bit cleaner than that of PRCE method (see Fig. 6).    Observing the stability diagram of the SSI-cov method (see Fig. 9), one can easily obtain the modal frequencies of 0.91 Hz, 3.06 Hz and 22.57 Hz. In contrast, no frequencies can be firmly determined from the stability diagram of the ERA method (see Fig. 8). Though the poles associated with the frequencies of 0.91 Hz and 3.06 Hz are struggling to spring out from the model order of 6, they do not get stable even when the model order gets 40.
In summary, the comparison study between a deterministic modal identification method and its counterpart in stochastic modal analysis shows that the stochastic methods are more suitable for a real platform that is excited by natural loading, such as the floating ice. By observing the stability diagrams from the stochastic methods, more structural frequencies can be estimated. Besides, the poles can get stable at a lower model order. Their stability diagrams are also cleaner than the stability diagrams from the deterministic methods. The reason is that for field experiment, there exists much noise from numerous sources. Only in stochastic modal analysis methods, noises are considered in their models and handled in their algorithms.

High-order methods versus low-order methods
As mentioned, all time domain methods need to adopt a mathematical model, either a high-order matrix polynomial differential/difference equation or a first-order state-space equation, to characterize the system's dynamic behavior. Mathematically, any high-order polynomial equation can be easily converted into a first order state-space form. However, the reverse is not always true. In this section, the comparison study is between high-order matrix polynomial and first-order state-space model methods. For the deterministic modal identification methods, the PRCE method (based on a high-order matrix polynomial model) and the ERA method (based on a first-order state-space matrix model) are compared. For the stochastic modal identification methods, the IV method (based on a high-order ARMA model) and the SSI-cov method (based on a stochastic state-space model) are compared.
For both the PRCE and ERA methods, it is hard to extract the modal frequencies while the stability diagrams are constructed with the model order ranging from 2 to 40 (see Figs. 6 and 8). Thus, making a comparison based on Figs. 6 and 8 is questionable. For a better illustration, to be compared are the stability diagrams with a model order ranging from 2 to 50 (see Figs. 10 and 11). By observing these figures, the PRCE method can roughly estimate the frequencies of 3.06 Hz and 22.57 Hz, while the ERA method can roughly identify the frequencies of 0.91 Hz and 3.06 Hz. For the mode which can be roughly identified by both methods, i.e. the frequency of 3.06 Hz, it is obvious that implementing the PRCE method requires a higher model order to form stable poles.
Through observing Fig. 7 and Fig. 9, the frequencies of 0.91 Hz, 3.06 Hz and 22.57 Hz can be identified by the SSIcov method, whereas the IV method can only identify the frequencies of 3.06 Hz and 22.57 Hz. Besides, for the identifiable frequencies, the poles from the SSI-cov method can become stable at a lower model order.
The comparison study between the high-order and loworder methods demonstrates that the EMA methods based on first-order state-space models can identify more modal information from the response of a real platform excited by floating ices. The superiority of the first-order methods can be explained by the numerical conditioning. It has been theoretically illustrated that methods based on first-order statespace models are more likely to be well-conditioned (with a smaller conditioning number) (Hu et al., 2014), and thus can generate better numerical performance than those based on high-order polynomial models.

SSI-cov method versus SSI-data method
The foregoing study has demonstrated that the stochastic methods are superior to the deterministic methods for the identification problem in this paper. Furthermore, the firstorder methods outperform the high-order methods. Thus, a comparison study between the two first-order stochastic methods, the SSI-cov method and the SSI-data method, will  YANG Wen-long, WANG Shu-qing China Ocean Eng., 2020, Vol. 34, No. 4, P. 558-570 567 be attractive.
The SSI-cov method and the SSI-data method are similar in many aspects. Firstly, they are modeled by the same state-space equations, and utilize exactly the same noisy response data. Secondly, they both estimate the modal parameters through a system realization algorithm that involves the truncated singular value decomposition of a kernel matrix: (1) the projection of the row space of the future outputs into the row space of the past outputs for the SSI-data method, and (2) the output covariance Toeplitz matrix for the SSI-cov method. After the corresponding matrix being formed or computed, the implementation of these methods is similar. It consists of computing the SVD of the kernel matrix, truncating the SVD to the model order n, estimating the observability matrix and controllability matrix (or Kalman filter state space for SSI) by splitting the SVD in two parts, and finally estimating system matrices from the observability matrix. In practical applications, the sole difference of the two methods is related to the computation of either projection matrix (see Eq. (28)) or Toeplitz matrix (see Eq. (23)) from the same noisy raw data, noting that both the projection and covariance operations are aimed to cancel out the (uncorrelated) noise. The projection employed in the SSI-data method and the output covariance in the SSI-cov method are also closely related to each other theoretically. They are in similarity transformation when the number of rows i is the same.
While choosing i = 256, the stability diagram from the SSI-cov method is shown in Fig. 12, which is strikingly different from Fig. 5, in particular: (1) The stability diagram from the SSI-data method is much cleaner than that from the SSI-cov method. (2) The weakly excited pole near 7.67 Hz exists in the diagram of SSI-data, but none in that of SSIcov. (3) While only one strongly excited pole near 3.06 Hz and 22.57 Hz is shown in the diagram of SSI-data method, multiple close poles are stabilized in that of SSI-cov.
From the comparison of the stabilization diagrams, there is no doubt that the SSI-data method is superior to the SSIcov method. The SSI-cov and SSI-data methods are based on the same stochastic model and developed by using system realization theory. The main difference between these methods is the kernel matrix employed to obtain the observability matrix associated with the system. The numerical performances of these two methods must originate from the numerical algorithm associated with those kernel matrices. Implementing the SSI-cov method is likely to mistakenly identify spurious (computational) poles to be true (physical) poles from using stability diagrams. This important finding is from the observation that the spurious poles could be stabilized in its stability diagram at frequency locations that are near the poles associated with strongly excited modes. The stabilization of spurious poles in the covariance-driven methods can be attributed to the "systematic error" caused by the round-off error in the (square operation) computation of the correlation functions or covariance matrices. No spurious poles would be stabilized when the numerical implementation of the SSI-data method has computed the required projection through a robust QR decomposition algorithm. Furthermore, the SSI-data method also has a better capability to identify the weakly excited modes. A possible explanation for why the SSI-cov method failing to identify the weakly excited mode of 7.67 Hz is that the SSI-cov method squares up the original signal and the square operation on a small number will become even a smaller number.

Concluding Remarks
In this paper, the possibility of employing the response from natural ice loading for identifying modal parameters of offshore platforms is substantiated. A field measurement on a real jacket platform was conducted and the acceleration data excited by floating ices were measured and utilized. From two segments, the modal frequencies and damping ratios of four structural modes can be successfully extracted. The estimated frequencies from both segments are almost identical, which demonstrates that the identification results are trustworthy.
Furthermore, by taking the test platform as a benchmark object, the numerical performance of five popular time-domain modal analysis methods was systematically compared. With the stability diagrams, it was found that those methods have different capabilities for identifying modal parameters: (1) The comparison on the PRCE method versus the IV method and the ERA method versus the SSI-cov method demonstrated that the stochastic methods have much better performance than their counterparts in deterministic modal analysis. The stability diagrams from the stochastic methods were much cleaner than those from the deterministic methods and more structural frequencies can be easily estimated. One can conclude that the stochastic methods are more suitable for a real platform that is excited by natural loading, such as the floating ice. The reason is that in the stochastic modal analysis methods, noises are considered in their models and handled in their algorithms.
(2) Compared with the high-order methods (the PRCE method and the IV method), the first-order methods (the ERA method, the SSI-cov method and the SSI-data method) have better numerical performance. More modal parameters can be extracted from the first-order methods. The poles associated with the structural frequencies from the first-order methods can get more stable at a lower model order than that of the high-order methods.
(3) By the measured data from a real offshore platform, one can draw the conclusion that the SSI-data method outperforms the SSI-cov method when they were employed together with a stability diagram to identify the true system poles. Implementing the SSI-cov method was likely to form spurious stable poles in its stability diagram at the frequency locations close to poles associated with strong modes. Besides, the SSI-data method also had a better capability to identify the weakly excited modes.