

\section{Introduction}

As was detailed in Chapter~\ref{chap:cr}, a key function for a practical assessment of the
extreme value distribution of stochastic response processes, is the average rate of upcrossings of high levels by the response. An important class of such response processes can be expressed as a second order stochastic Volterra series\index{Volterra series}\index{stochastic Volterra series}, that is, a stochastic Volterra series that has been truncated after the second order term \citep{sche80}. A substantial amount of work has been done to derive methods for efficient analysis of this model, starting with the seminal
paper by \cite{kasi47}. Later, with the development of the offshore industry, this paper had an impact on investigations of the response statistics of large floating structures. Early contributions in this field of research were made, among many others, by \cite{neal74}, \cite{vinj83}, \cite{lang84}, \cite{naes85jsr,naes90pem} and \cite{dosp90}.

The type of stochastic Volterra series models that will be studied
here, can be expressed as a sum of a linear and a nonlinear,
quadratic transformation of a Gaussian process. Such a
representation of the response process would apply to the standard
model for expressing the total wave forces or horizontal excursion
responses of e.g. a tension leg platform in a random sea way. It
also applies to the response of a linear structure to a quadratic
wind loading where the wind speed is modeled as a Gaussian
process. It would also apply to the representation of large
deformations of flexible structures subjected to Gaussian loads in
which the strains/stresses exhibit a strong quadratic effect. The
problem of determining the marginal probability distribution
function of such response processes has been solved in the sense
that computer programs are available that allow very accurate
numerical calculations of these functions, cf. 
\cite{najo92jomae}. 

The focus of this chapter is to show that it is also possible to accurately calculate the
average upcrossing rate of second order stochastic Volterra models. This allows the formulation of approximate extreme value distributions for such response processes. An interesting aspect of the development in this chapter is the surprising structural complexity of the problem of calculating the upcrossing rate of a second order Volterra model. It is also hoped that this chapter may serve to illustrate the power of the point process approach to practical extreme value analysis via the upcrossing rate function.  The performance of the numerical method that has been developed will be illustrated by application to three specific examples.


\section{The Response Process}\index{response process}

As already stated, the object of study in this chapter is a
stochastic response process $Z(t)$ modeled as a second order
stochastic Volterra series. Specifically, it is assumed that
$Z(t)$ can be written as a sum of a linear, first order response
component $Z_1(t)$ and a nonlinear, second order component
$Z_2(t)$. That is, cf. \cite{naes90pem},
\begin{equation}\label{eq:cf1}
Z(t) = Z_1(t) + Z_2(t),
\end{equation}
where
\begin{equation}\label{eq:cf2}
Z_1(t) = \int_0^\infty h_1(\tau) X(t-\tau) \, d \tau,
\end{equation}
and
\begin{equation}\label{eq:cf3}
Z_2(t) = \int_0^\infty \int_0^\infty h_2(\tau_1, \tau_2)
X(t-\tau_1) X(t-\tau_2) \, d \tau_1 d \tau_2.
\end{equation}

In Eqs.~\eqref{eq:cf2} and \eqref{eq:cf3}, $X(t)$ denotes a stationary, real
Gaussian process. $X(t)$ could represent a random wave elevation
process or a stochastic wind velocity field. It was chosen here
to limit the exposition to the case of a unidirectional situation.
How to deal with the multi-directional case, is explained in detail
in \cite{naes90pem}. Based on the results from this reference, it
will be recognized that all results obtained in this chapter apply
equally well to the multi-directional case.

The functions $h_1(\tau)$ and $h_2( \tau_1, \tau_2)$ characterize
the physical system that is modeled. $h_1(\tau)$ is an ordinary
impulse response function defining a linear dynamical system.
$h_2( \tau_1, \tau_2)$, which is referred to as the quadratic
impulse response function, characterizes the second-order
properties of the physical system, but in contrast to the linear
impulse response function, it does not have a direct physical
interpretation.

To derive expressions suitable for practical numerical
calculations, the input process $X(t)$ is represented as follows,
\begin{equation}\label{eq:cf4}
X(t) = \sum_{j=-N}^N \Big[S_X(\omega_j) \Delta \omega \Big]^{1/2}
B_j \, e^{\ii \omega_j t},
\end{equation}
where $S_X(\omega)$ denotes the two-sided spectral density of
$X(t)$. Throughout this chapter, when the summation index runs from
negative to positive values, it invariably omits zero. $ 0 <
\omega_1 < \cdots < \omega_N $ is an equidistant discretization of
the pertinent part of the positive frequency axis. $\omega_{-i}= -
\omega_i$, $\Delta \omega = \omega_{i+1} - \omega_i$. The
assumption of an equidistant discretization is adopted for
simplicity of presentation and is not necessary. In fact, often a
non-equidistant version is used to avoid having too many
frequencies, which is sometimes convenient. The formulas are
easily adapted to cover the situation of non-equidistant
discretization. $ \{ B_{i} \}$ is a set of independent, complex
Gaussian $N(0,1)$-variables with independent, identically
distributed real and imaginary parts. These variables can be
assumed to satisfy the relation $ B_{-i} = B_{i}^*$, where
* signifies complex conjugation. $\ii^2 = -1$. Hence, $X(t)$ becomes a Gaussian process with
zero mean value. The case of a non-zero mean value does not create
any difficulties for the analysis, as is easily recognized.

By substituting Eq.~\eqref{eq:cf4} into Eqs.~\eqref{eq:cf2} and \eqref{eq:cf3}, the
following expressions are obtained,
\begin{equation}\label{eq:cf5}
Z_1(t) = \sum_{i=-N}^{N} q_{i} B_{i}
e^{\ii \omega_i t},
\end{equation}
where\label{eq:cf6}
\begin{equation}
q_{i} = \hat{H}_1(\omega_i)[S_X(\omega_i) \Delta \omega ]^{1/2},
\end{equation}
and
\begin{equation}\label{eq:cf7}
Z_2(t) = \sum_{i=-N}^{N} \sum_{j=-N}^{N} Q_{ij} B_{i} B_{j}^*
e^{\ii (\omega_i - \omega_j)t},
\end{equation}
where
\begin{equation}\label{eq:cf8}  
Q_{ij} = \hat{H}_2(\omega_i, -\omega_j)
\cdot [S_X(\omega_i) S_X(\omega_j)]^{1/2} \Delta \omega.
\end{equation}

The function $\hat{H}_1(\omega)$ denotes the linear transfer
function\index{linear transfer function}, corresponding to $h_1(\tau)$. Typically, the second
order Volterra series is used for modeling the stochastic loading
process. For the response process to be of the same type, it is
necessary that the equations of motion of the dynamic system
considered are linear and time invariant. Such a system is
characterized by a linear transfer function, $\hat{L}(\omega)$
say. The implication of this is that $\hat{H}_1(\omega)$ can be
expressed as follows,
\begin{equation}\label{eq:cf9}
\hat{H}_1(\omega) = \hat{L}(\omega) \hat{K}_1(\omega).
\end{equation}
$\hat{L}(\omega)$ is assumed to be given as follows,
\begin{equation}\label{eq:cf10}
\hat{L}(\omega ) = (-\omega^2 M + \ii \omega C + K)^{-1},
\end{equation}
where $M=M(\omega)$, $C=C(\omega)$  and $K$ are appropriate mass,
damping and stiffness parameters. The function $\hat{K}_1(\omega)$
is a linear transfer function associated with the loading process.

The function $\hat{H}_2(\omega,\omega')$ denotes the quadratic
transfer function (QTF)\index{quadratic transfer function}, which depends on two frequencies. It is
obtained as a Fourier transform of $h_2(\tau_1, \tau_2)$.
$\hat{H}_2(\omega,\omega')$ can be expressed as follows,
\begin{equation}\label{eq:cf11}
\hat{H}_2(\omega,\omega') =
\hat{L}(\omega + \omega' ) \hat{K}_2(\omega,\omega').
\end{equation}
The function $\hat{K}_2(\omega,\omega')$  denotes a QTF for the
quadratic forces on the structure.

It is shown by \cite{naes87jsv}, see also 
\cite{joko70}, that by solving the eigenvalue problem (assumed
nonsingular),
\begin{equation}\label{eq:cf12}
Q v_j = \lambda_j v_j,
\end{equation}
to find the eigenvalues $\lambda_j$ and orthonormal eigenvectors
$v_j$, $j= - N, \ldots, -1,$ $1, \ldots, N$, the quadratic
response process can be represented as
\begin{equation}\label{eq:cf13}
Z_2(t) = \sum_{j=-N}^N \lambda_j \, W_j(t)^2.
\end{equation}
Here $W_j(t)$, $j=-N, \ldots, -1,$ $1, \ldots, N$ are real
stationary Gaussian $N(0,1)$-processes which can be represented as
follows
\begin{equation}\label{eq:cf14}
W_j(t) = \sum_{k=-N}^N v_j(\omega_k) B_k e^{\ii \omega_k t},
\end{equation}
where $v_j(\omega_k)$ denotes the $k$th component of $v_j$. Note
that the property $v_j(\omega_{-k})=v_j(\omega_k)^*$ can be
assumed, cf. \cite{naes90pem}. For each fixed $t$, $\{W_j(t)\}$
becomes a set of independent Gaussian variables.

Having achieved the desired representation of the quadratic
response $Z_2(t)$, it can then be shown that the first order
response can be expressed as
\begin{equation}\label{eq:cf15}
Z_1(t) =  \sum_{j=-N}^{N} \beta_j W_j(t).
\end{equation}
The (real) parameters $\beta_j$ are given by the relations
\begin{equation}\label{eq:cf16}
\beta_j = \sum_{k=-N}^N \hat{H}_1^*(\omega_k) \sqrt{S_X(\omega_k)
\Delta \omega}\, v_j(\omega_k).
\end{equation}

This gives us the following representation of the total response
process \citep{vinj83,naes85jsr}
\begin{equation}\label{eq:cf17}
Z(t) =  \sum_{j=-N}^{N} \{\lambda_j \, W_j(t)^2 + \beta_j W_j(t)
\}.
\end{equation}

Based on this representation, \cite{naes87jsv} describes how to
calculate the statistical moments of the response process.
Asymptotic upper and lower bounds on the mean level upcrossing
rate are also given, but these bounds are valid only if specific
conditions on the ratios $\beta_j/\lambda_j$ are satisfied. It
turns out that these conditions are very often violated in
practical applications. Our goal is therefore to find a numerical
method for accurate calculation of the mean level upcrossing rate
of the total response process $Z(t)$. This is the topic of the
next section.


\section{The Average Crossing Rate }

The calculation of the mean upcrossing rate\index{mean upcrossing rate} $\nu_Z^+(\cdot; t)$ at time $t$ of a
stochastic process $Z(t)$ is usually based on the Rice formula. It
states that
\begin{equation}\label{eq:cf18}
\nu_Z^+(\zeta; t) = \int_0^\infty \dot{z} f_{Z(t)\dot{Z}(t)}(\zeta,
\dot{z}) \, d \dot{z},
\end{equation}
where $f_{Z(t)\dot{Z}(t)}(\cdot, \cdot)$ denotes the joint probability
density function of $Z(t)$ and $\dot{Z}(t) = dZ(t)/dt$.  A
direct application of Eq.~\eqref{eq:cf18} requires the calculation of
$f_{Z(t)\dot{Z}(t)}$. An alternative approach, which turns out to be
ideally suited for the problem at hand, is to express the mean
crossing rate in terms of the characteristic function
$M_{Z\dot{Z}}(u, v) = E[\exp\{\ii (u\, Z(t) + v \, \dot{Z}(t) )\} ]$. It
has been shown by \cite{nar200} that in general,
\begin{eqnarray} \label{eq:cf19}
 \nu_Z^{+}(\zeta, t) = &-& \frac{1}{(2 \pi)^2}
        \int_{-\infty}^{+\infty}
        \int_{-\infty}^{+\infty}
        \frac{1}{v}
        \frac{\partial M_{Z\dot{Z}}(u, v)}{\partial v}\,
        e^{-\ii \zeta \, u}\, d u
        \, d v \nonumber \\
         &-& \frac{\ii}{4 \pi}
        \int_{-\infty}^{+\infty}
        \left(\frac{\partial M_{Z\dot{Z}}( u, v)}{\partial v}\right)_{v=0} \,
        e^{-\ii \zeta \,u}\, d u,
\end{eqnarray}
where the integral wrt $v$ is interpreted as a principal value
integral in the following sense:
$\int_{-\infty}^\infty=\lim_{\varepsilon \to 0}(
\int_{-\infty}^{-\varepsilon}+\int_{\varepsilon}^\infty)$. A
heuristic derivation of this formula is given in Section~\ref{Ap1}.

For the stationary case, it can be shown that the last integral on
the right hand side of Eq.~\eqref{eq:cf19} vanishes. Hence, for a
stationary process, the following formula is obtained
\citep{vinj83,nar200},

\begin{equation}\label{eq:cf20}
\nu_Z^{+}(\zeta) = - \frac{1}{(2 \pi)^2}
        \int_{-\infty}^{+\infty}
        \int_{-\infty}^{+\infty}
        \frac{1}{v}
        \frac{\partial M_{Z\dot{Z}}(u, v)}{\partial
        v}\,
        e^{-\ii \zeta \,u}\, d u
        \, d v,
\end{equation}
where the upcrossing rate is now independent of $t$, and again, the integral with respect to $v$ is to be interpreted as a principal value integral as described above.

An alternative expression useful for numerical calculations of the upcrossing rate $\nu_Z^+(\cdot)$, whether the process is stationary or not, can be obtained by considering the characteristic 
function as a function of two complex variables. It can then often
be shown that this new function becomes holomorphic in suitable
regions of $\mathbf{C}^2$, where $\mathbf{C}$ denotes the complex
plane. Under suitable conditions, the use of complex function theory allows the
derivation of the following alternative expression for the crossing rate, cf. \cite{naka04pem}.
\begin{align}\label{mcr7}
\nu_Z^+(\zeta)= -\,\frac{1}{(2\, \pi)^2}  \int_{-\, \infty - \ii a
}^{\infty - \ii a} \int_{- \infty - \ii b}^{\infty - \ii b}
\frac{1}{w^2}\, M_{Z\dot{Z}}( z, w)\, e^{-\, \ii z \zeta} \, d z \, d w,
\end{align}
where $0 < a < a_1$ for some positive constant $a_1$, and $b_0 < b
< b_1$ for some constants $b_0 < 0$ and $b_1 > 0$. 



The calculation of the characteristic function\index{characteristic function} $M_{Z\dot{Z}}$ is
discussed in \cite{grig95}, but no explicit expression is derived.
Here, a different approach will be followed, which leads to a
convenient explicit representation of the characteristic function
suitable for calculation of the integrals appearing in Eqs.~\eqref{eq:cf19}, \eqref{eq:cf20} or \eqref{mcr7}. 
To this end, consider the multidimensional Gaussian
vectors $W=(W_{-N}, \ldots, W_{N})'$ and $\dot{W}=(\dot{W}_{-N},
\ldots, \dot{W}_{N})'$. It is obtained that the covariance matrix
of $(W', \dot{W}')'$ is given by,
\begin{equation}\label{eq:cf21}
 \Sigma= \left( \begin{array}{cc} \Sigma_{11} & \Sigma_{12} \\
         \Sigma_{21} & \Sigma_{22} \end{array}\right),
\end{equation}
where $\Sigma_{11}= I$ = the $2N\times2N$ identity matrix,
$\Sigma_{12}=(r_{ij})= (E[W_i\dot{W}_j])$,
$\Sigma_{21}=(E[\dot{W}_i W_j])$ and $\Sigma_{22}=(s_{ij})=
(E[\dot{W}_i\dot{W}_j])$; $i, j=-N, \ldots, -1,$$1, \ldots, N$.
$r_{ij}= - \,r_{ji}$ and $\Sigma_{12}= \Sigma_{21}'$. It follows
from Eq.~\eqref{eq:cf14}, that the entries of the covariance matrix
$\Sigma$ can be expressed in terms of the eigenvectors $v_j$, cf.
\cite{naes87jsv}. Let $\Lambda=\diag(\lambda_{-N}, \ldots,
\lambda_{N} )$ be the diagonal matrix with the parameters
$\lambda_j$ on the diagonal, and let $\beta=(\beta_{-N}, \ldots,
\beta_{N})'$, cf. Eq.~\eqref{eq:cf16}. It will be shown in Section~\ref{Ap2} that \citep{nar100},
\begin{equation}\label{eq:cf22}
 M_{Z\dot{Z}}(u, v)= \frac{1}{\sqrt{\det(A)}}\,\, \exp \Bigg\{-\frac{1}{2}\, v^2\,
                       \beta'\,C \,\beta\,+
                   \frac{1}{2}\,d'\,A^{-1}\,d
               \Bigg\},
\end{equation}
where
\begin{equation}\label{eq:cf23}
A = A(u, v)= I-2 \ii \, u \,\Lambda\, -2\,\ii \,v\,\big(\Lambda\,
               \Sigma_{21}+\Sigma_{12}\,\Lambda \big) + 4\,v^2 \Lambda\, C \,\Lambda,
\end{equation}
\begin{equation}\label{eq:cf24}
C = \Sigma_{22}-\,\Sigma_{21}\, \Sigma_{12},
\end{equation}
\begin{equation}\label{eq:cf25}
d=d(u, v) = \Big(\ii \,u\,I + \ii \,v\,\Sigma_{12}- 2\,v^2\,
\Lambda\,C\,\Big)\beta. 
\end{equation}
 


\section{Numerical Calculation}

Early efforts to carry out numerical calculation of the mean
crossing rate using Eq.~(\ref{mcr7}) has been reported in
\cite{naka04pem}. These initial investigations indicated that the
method had the potential to provide very accurate numerical
results. Eq.~(\ref{mcr7}) will be rewritten as follows,
\begin{flalign}\label{nc0}
&\nu_Z^+(\zeta)= -\,\frac{1}{(2\, \pi)^2}  \int_{-\, \infty - \ii
	a }^{\infty - \ii a} \frac{1}{w^2}\, I( \zeta, w) \, d w,
\end{flalign}
where
\begin{flalign}\label{nc1}
I&=I(\zeta, w) = \int_{- \infty - \ii b}^{\infty - \ii b} \, M(
z, w)\, e^{-\, \ii z \zeta} \, d z \nonumber\\
&=\int_{- \infty - \ii b}^{\infty - \ii b} \exp \{ -\, \ii z \zeta
+ \ln M( z, w) \} \, d z.
\end{flalign}


A numerical calculation of the mean upcrossing rate can start by
calculating the function $I(\zeta, w)$ for specified values of
$\zeta$ and $w$. However, a direct numerical integration of
Eq.~(\ref{nc1}) is made difficult by the oscillatory term
$\exp\{ - \ii \Re(z) \zeta \}$, where $\Re(z)$ denotes the real part of the complex number $z$. This problem can be avoided by
invoking the method of steepest descent, also called the saddle
point method. For this purpose, write,
\begin{flalign}\label{nc2}
g(z)&= g(z; w) = -\, \ii z \zeta + \ln M( z, w) \nonumber\\
&=\phi(x,y) + \ii \, \psi(x,y),
\end{flalign}
where $z = x + \ii y$. $\phi(x,y)$ and $\psi(x,y)$ become real
harmonic functions when $g(z)$ is holomorphic. The idea is to
identify the saddle point of the surface $(x,y) \to \phi(x,y)$
closest to the integration line from $- \infty - \ii b$ to $\infty
- \ii b$. By shifting this integration line to a new integration
contour that passes through the saddle point, and then follows the
path of steepest descent away from the saddle point, it can be
shown that the function $\psi(x,y)$ stays constant, and therefore
the oscillatory term in the integral degenerates to a constant.
This is a main advantage of the method of steepest descent for
numerical calculations. It can be shown that the integral does not
change its value as long as the function $g(z)$ is a holomorphic
function in the region bounded by the two integration contours and
if the integrals vanish along the contour segments required to
close the region.

If $z_s$ denotes the identified saddle point, where $g'(z_s)=0$,
the steepest descent path away from the saddle point will follow
the direction given by $- \, {g'(z)}^*$, for $z \ne z_s$, cf.
\cite{henr77}. Typically, the singular points of the function $g$
will be around the imaginary axis, which indicates that the
direction of the paths of steepest descent emanating from the
saddle point will typically not deviate substantially from a
direction orthogonal to the imaginary axis. This provides a guide
for setting up a numerical integration procedure based on the path
of steepest descent. First the saddle point $z_s$ is identified.
Then the path of steepest descent starting at $z_s$ and going
'right', is approximated by the sequence of points
$\{z_j\}_{j=0}^\infty$ calculated as follows:
\begin{equation}\label{nc3}
z_0 = z_s, \qquad z_1 = z_s + h,
\end{equation}
\begin{flalign}\label{nc4}
\Delta z_j = - \, \frac{{g'(z_j)}^*}{|g'(z_j)|}\,h, \quad j=1,2,
\ldots
\end{flalign}
\begin{equation}\label{nc5}
z_{j+1} = z_j + \Delta z_j, \quad j=1,2, \ldots
\end{equation}
where $h$ is a small positive constant.

Similarly, the path of steepest descent going 'left' is
approximated by the sequence $\{z_j\}_{j=0}^{-\infty}$ calculated
by,
\begin{equation}\label{nc7}
z_{-1} = z_s - h,
\end{equation}
\begin{flalign}\label{nc8}
\Delta z_j = - \, \frac{{g'(z_j)}^*}{|g'(z_j)|}\,h, \quad j=-1,-2,
\ldots
\end{flalign}
\begin{equation}\label{nc9}
z_{j-1} = z_j + \Delta z_j, \quad j=-1,-2, \ldots
\end{equation}

A numerical estimate $\hat{I}$ of $I$ can be obtained as follows,
\begin{equation}\label{nc10}
\hat{I} = \hat{I}^+ + \hat{I}^-,
\end{equation}
where
\begin{flalign}\label{nc11}
\hat{I}^+ = \frac{h }{2}\,\exp \{g(z_s)\}  +  \sum_{ j=1}^K \Delta
z_j \, \exp \{ g(z_j)\},
\end{flalign}
and
\begin{flalign}\label{nc12}
\hat{I}^- = \frac{h }{2}\,\exp \{g(z_s)\}  -  \sum_{ j=-1}^{-K}
\Delta z_{j} \, \exp \{ g(z_{j})\},
\end{flalign}
for a suitably large integer $K$.

A numerical estimate $\hat{\nu}_Z^+(\zeta)$ of the mean crossing
rate can now be calculated by the sum,
\begin{flalign}\label{nc13}
\hat{\nu}_Z^+(\zeta)= -\,\frac{1}{(2\, \pi)^2}\, \Re \Big\{ \sum_{
	j=-L}^{L} \frac{1}{w_j^2} \, \hat{I}(\zeta, w_j ) \, \Delta w_j
\Big\},
\end{flalign}
where the discretization points $w_j$ are chosen to follow the
negative real axis from a suitably large negative number up to a
point at $-\varepsilon$, where $0 < \varepsilon \le a$, then
follow a semi-circle in the lower half plane to $\varepsilon$ on
the positive real axis, and finally follow this axis to a suitably
large positive number. Since the numerical estimate does not
necessarily have an imaginary part that is exactly equal to zero,
the real part operator $\Re$ has been applied.

Generally, the CPU time required to carry out the computations
above can be quite long, depending on the size of the problem,
which is related to the number $N$ of eigenvalues. It is therefore
of interest to see if approximating formulas are accurate enough.
The first such approximation to have a look at, is the
Laplace approximation for the inner integral over the saddle point
\citep{henr77}. The simplest version of this approximation, adapted
to the situation at hand, leads to the result,
\begin{flalign}\label{nc14}
I=I(\zeta, w) \approx \sqrt{ \frac{2 \pi}{ - \frac{\partial^2
			g(z_s;w)}{\partial x^2} } } \, \exp\{\,g(z_s;w)\},
\end{flalign}
which can be substituted directly into Eq.~(\ref{nc13}), leading
to an approximation of ${\nu}_Z^+(\zeta)$, which is denoted by
$\tilde{\nu}_Z^+(\zeta)$.

This approximation can also be exploited in the following way: 1) The
full method is used for an inner interval of $w$-values, which
contribute significantly to the integral in Eq.~(\ref{nc0}). 2) A
Laplace approximation is then used in an outer interval of
$w$-values where the contribution is less than significant. Of
course, the level of significance is chosen according to some
suitable criterion. By this procedure, the CPU time was reduced by
factor of about 3. This method will be referred to as the hybrid
method, and the corresponding approximation of ${\nu}_Z^+(\zeta)$
is denoted by $\check{\nu}_Z^+(\zeta)$.

A simple approximation proposed in \cite{tena99,tena03} is worth a
closer scrutiny. It is based on the widely adopted simplifying
assumption that the displacement process is independent of the
velocity process. This leads to an alternative approximation of
${\nu}_Z^+(\zeta)$, which is denoted by
$\overline{\nu}_Z^+(\zeta)$. It is given by the formula,
\begin{flalign}\label{nc15}
\overline{\nu}_Z^+(\zeta)= \nu_Z^+(\zeta_{ref})
\frac{f_Z(\zeta)}{f_Z(\zeta_{ref})},
\end{flalign}
where $f_Z$ denotes the marginal probability density of the surge response, and
$\zeta_{ref}$ denotes a suitable reference level, typically the
mean response. Here, $\zeta_{ref}$ has been chosen as the point
where $f_Z$ assumes its maximum, which corresponds well with the
mean response level. A general approximation for
$\nu_Z^+(\zeta_{ref})$ is given in \cite{tena03}. If only
slow-drift response is considered, a good approximation is
obtained by putting $\nu_Z^+(\zeta_{ref}) \approx 1 / T_0$, where
$T_0 = 2 \pi / \omega_0$ is the slow-drift period. The advantage
of Eq.~(\ref{nc15}) is that the rhs is much faster to calculate
than the exact formula.

A few comments on why the mean upcrossing rate has practical significance may seem appropriate. As is already known, the extreme value distribution is not completely determined by the mean upcrossing rate. This is true only when the upcrossing events of the high response levels can be
assumed to be statistically independent. Usually that is a good
approximation except when the total damping is very small. For such
cases, \cite{naes99jomae} has developed a simple, but effective, method
to account for the effect of low damping on the extreme value
distribution. This method is based on the mean upcrossing rate and the appropriate damping parameter. This is often advantageous since no time series of response is required. However, when those are available, the key to an accurate estimation of the extreme value distribution would be the application of the ACER method. 





\section{Numerical Examples}

\subsection{Slow-drift response}\index{slow-drift response}

To illustrate the accuracy of the numerical method, the first
example concerns a simple model for the slow-drift response of a
moored offshore structure, cf. \cite{nama00}. Specifically, the
response process for this case may be written as
($\lambda=\lambda_1=\lambda_2$),
\begin{equation}\label{eq:cf31}
Z_2(t)=\lambda \bigg( W_1(t)^2 + W_2(t)^2 \bigg).
\end{equation}
For this process, the upcrossing rate $\nu_{Z_2}^+(z)$ is given by
the relation,
\begin{equation}\label{eq:cf32}
\nu_{Z_2}^+(\zeta)=\frac{\hat{\sigma}_1}{\sqrt{2 \pi }} \exp \Big( -
\frac{\zeta}{2 \lambda} + \frac{1}{2} \ln \big(\frac{\zeta}{\lambda}\big)
\Big),
\end{equation}
where $\hat{\sigma}_1 = \sqrt{s_{11}-(r_{12})^2}$. This special
case provides a suitable test for the accuracy of the numerical
method.

Let $\tilde{\nu}_{Z_2}^+(\zeta)$ denote the mean upcrossing rate of
$Z_2(t)$ calculated by the numerical method. Table~\ref{tab:7.1} compares the
analytical with the numerical upcrossing rate for different
levels, and it is seen that the agreement is very good indeed.


\begin{table}[!htb]
	\begin{center}
		\caption{Comparison of analytical and numerical mean upcrossing rate. \label{tab:7.1}}
		\vspace{2mm}
		\begin{tabular}{c c c}
			\hline 
			$\zeta$   &    $\nu_{Z_2}^+(\zeta)$   &  $\tilde{\nu}_{Z_2}^+(\zeta)$\\ 
			\hline  
			0.5 & $1.37 \cdot 10^{-2}$ & $1.37 \cdot 10^{-2}$\\ 
			 
			1.0 & $4.31 \cdot 10^{-3}$ & $4.30 \cdot 10^{-3}$ \\ 
			 
			2.0 & $2.995 \cdot 10^{-4}$ & $2.995 \cdot 10^{-4}$\\ 
			 
			3.0 & $1.802 \cdot 10^{-5}$ & $1.793 \cdot 10^{-5}$\\ 
			 
			5.0 & $5.618 \cdot 10^{-8}$ & $5.618 \cdot 10^{-8}$\\ 
		
			7.0 & $1.605 \cdot 10^{-10}$ & $1.592 \cdot 10^{-10}$\\ 
			\hline
		\end{tabular}
	\end{center}
\end{table}




\subsection{Moored deep floater}\label{mdf}\index{moored deep floater}

The numerical results presented in this example are based on a
specific model structure \citep{nakate06aor}. It is a moored deep floater (MDF), which is also called a spar-buoy, with main particulars as detailed in Table~\ref{tab:mdf}.
Figure~\ref{fig:mdf} shows the submerged part of the floater in
the form of a computer mesh, which is used for the calculation of
the hydrodynamic transfer functions. The total mass (including
added mass) of the MDF is $M = 12.5 \cdot 10^6$ kg. The damping
ratio is set equal to $\xi=$ 0.06, and the natural frequency in
surge or sway is $\omega_0$ = 0.047 rad/s.  Note that the second order
theory is based on the assumption that the QTF
$\hat{H}_2(\omega_i, - \omega_j) = \hat{L}(\omega_i - \omega_j)\,
\hat{K}_2(\omega_i, - \omega_j)$, where $\hat{K}_2(\cdot, \cdot)$
is a QTF characterizing the slowly varying surge forces on the
MDF, and $\hat{L}(\cdot)$ is a linear transfer function for the
surge motion of the MDF, that is,
\begin{equation}\label{ex1}
\hat{L}(\omega)=\frac{1}{M \, [ - \omega^2 + 2 \ii \xi \omega_0
	\omega + \omega_0^2 ]}.
\end{equation}
An example of the quadratic transfer function for a floating offshore structure is presented in Subsection~\ref{ssec:tsrtlp}.


\begin{table}[!htb]
	\caption{Main particulars of the moored deep floater (MDF)}
	\begin{center}
		\label{tab:mdf}
		\begin{tabular}{l  r}
			% put some space after the caption
			& \\
			\hline
			Draught (m) & 80.0  \\
			Column diameter (m) & 10.0 \\
			Natural period surge/sway (s) & 133.5 \\
			Natural period yaw (s) & 121 \\
			\hline
		\end{tabular}
	\end{center}
\end{table}


\begin{figure}[!htb]
	\begin{center}
		\setlength{\unitlength}{0.0125in}%
		\dofig{5cm}{Figures/totsylc1.pdf}
	\end{center}
	\caption{Computer mesh of the submerged part of the moored deep
		floater and the near field of the sea surface.} \label{fig:mdf}
\end{figure}


The deep floater studied in this example could represent the supporting structure of a floating wind turbine. The Hywind turbine is a particular case of such a structure, where the concept of a moored deep floater is used as a supporting structure. A sketch of the Hywind turbine is presented in Fig.~\ref{fig:hywind}. 



\begin{figure}[!htb]
	\begin{center}
		\setlength{\unitlength}{0.0125in}%
		\dofig{5cm}{Figures/Hywind_illustration_g.pdf}
	\end{center}
	\caption{A sketch of the Hywind floating wind turbine (\copyright Equinor).} \label{fig:hywind}
\end{figure}


The random stationary sea state\index{stationary sea state} is specified by a JONSWAP
spectrum\index{JONSWAP spectrum}, which is given as follows,
\begin{equation}\label{mdfex2}
S_X(\omega)=\frac{\alpha g^2}{\omega^5} \exp\Bigg\{ - \frac{5}{4}
\Big(\frac{\omega_p}{\omega}\Big)^4 + \ln \gamma \, \exp\Big[ -
\frac{1}{2 \sigma^2} \Big(\frac{\omega}{\omega_p} - 1 \Big)^2
\Big] \Bigg\},
\end{equation}
where $g = 9.81$ ms$^{-2}$, $\omega_p$ denotes the peak frequency
in rad/s and $\alpha$,
$\gamma$ and $\sigma$ are parameters related to the spectral shape.
$\sigma = 0.07$ when $\omega \le \omega_p$, and $\sigma = 0.09$ when
$\omega > \omega_p$. The parameter $\gamma$ is chosen to be equal
to 3.0. The parameter $\alpha$ is determined from the following
empirical relationship,
\begin{equation}\label{mdfex3}
\alpha = 5.06 \Big(\frac{H_s}{T_p^2}\Big)^{2} \big( 1 - 0.287 \,
\ln \gamma \, \big).
\end{equation}
$H_s$ = significant wave height\index{significant wave height} and $T_p= 2 \pi/ \omega_p$ =
spectral peak wave period\index{spectral peak period}. For the subsequent calculations, $H_s$ = 10.0 m and $T_p$ = 12 sec. The natural frequency in surge is
0.047 rad/s, which is well below the range where the waves have
noticeable energy. This is why the second order, nonlinear term in
the Volterra expansion is needed to capture the resonant motions
in surge of the MDF.

To get an accurate representation of the response process, there
is a specific requirement that must be observed. Since the damping
ratio is only $6\%$, the resonance peak of the linear transfer
function for the dynamics is quite narrow. Hence, to capture the
dynamics correctly, the frequency resolution must secure a
sufficient number of frequency values over the resonance peak.
This usually leads to an eigenvalue problem with the $Q$-matrix of
size of the order of magnitude $100 \times 100$. Using the full
representation of this size in calculating the mean crossing rate
by the general method described here, would lead to very heavy
calculations. In order to reduce this, the effect of restricting the calculations by retaining only some of
the terms in Eq.~(\ref{eq:cf13}) has been investigated.

For the specific example considered, where exactly 100 (positive) frequencies have been used, the values of the obtained eigenvalues $\lambda_j$ have been plotted in
Figure~\ref{fig:eigenvalues}. It is seen that a substantial portion
of the response variance, which is given by $\Var[Z_2(t)] = 4 \,
\sum_{j=1}^N \lambda_j^2$, would be lost if only 10 or 20 eigenvalues
were retained. This is also a factor to consider when deciding on
the number of terms to retain.


\begin{figure}[!htb]
	\begin{center}
		\setlength{\unitlength}{0.012500in}%
		\dofig{6.5cm}{Figures/eig.pdf}
	\end{center}
	\caption{The 100 normalized eigenvalues $\lambda_j / \lambda_1$.}
	\label{fig:eigenvalues}
\end{figure}


In this example, the focus is on the slow-drift response. Hence,
only results for $Z_2(t)$ will be presented. In the tables,
$\hat{\nu}_{Z_2}^+(\zeta)$, $\check{\nu}_{Z_2}^+(\zeta)$,
$\tilde{\nu}_{Z_2}^+(\zeta)$, and $\overline{\nu}_{Z_2}^+(\zeta)$)
denote the mean upcrossing rate of $Z_2(t)$ calculated by the full
numerical method, the hybrid method, the Laplace approximation,
and the simplified method of Eq.~(\ref{nc15}), respectively.

To highlight the effect of the increment parameter $h$,
Table~\ref{tab:2} compares the results obtained by the full
numerical method for two values of $h$ for 10 eigenvalues, that
is, for a response representation retaining the first 10 terms.
The CPU time differs by a factor of roughly 10 between the two
choices of a value for $h$. Since the differences between the
calculated crossing rates are fairly small, the larger value was chosen to save CPU time.


\begin{table}[!htb]
	\caption{ Comparison of calculated mean upcrossing rate
		$\hat{\nu}_{Z_2}^+(\zeta)$ for different step lengths
		\label{tab:2}}
	\begin{center}
		\begin{tabular}{c c c}
			\hline
			$\eta = \zeta/ \lambda_1$   &    $h = 1.0 \cdot 10^{-3}$      & $ h = 1.0 \cdot 10^{-2} $\\
			\hline
			2.0 & $8.38 \cdot 10^{-3}$ & $ 8.38 \cdot 10^{-3}$\\
			5.0 & $3.93 \cdot 10^{-3}$ & $3.93 \cdot 10^{-3}$ \\
			10.0 & $5.53 \cdot 10^{-4}$ & $5.50 \cdot 10^{-4}$\\
			15.0 & $5.70 \cdot 10^{-5}$ & $5.65 \cdot 10^{-5}$\\
			20.0 & $5.34 \cdot 10^{-6}$ & $5.26 \cdot 10^{-6}$\\
			25.0 & $4.81 \cdot 10^{-7}$ & $4.71 \cdot 10^{-7}$\\
			\hline
		\end{tabular}
	\end{center}
\end{table}



Tables~\ref{tab:3} and \ref{tab:6} 
present the results obtained for 10 and 50
eigenvalues, respectively. It is apparent that there is some
variability of the calculated mean upcrossing rates depending on the
number of eigenvalues included in the analysis. Ideally, it would
therefore be desirable to carry out the calculations with at least
50 eigenvalues. 


\begin{table}[!htb]
	\caption{ Calculated mean upcrossing rates for 10 eigenvalues
		\label{tab:3}}
	\begin{center}
		\begin{tabular}{c c c c}
			\hline
			$\eta = \zeta/ \lambda_1$   &  $\hat{\nu}_{Z_2}^+(\zeta) $  & $\check{\nu}_{Z_2}^+(\zeta) $ & $\tilde{\nu}_{Z_2}^+(\zeta)$ \\
			\hline
			2.0 & $8.38 \cdot 10^{-3}$ & $ 8.38 \cdot 10^{-3}$ & $ 7.41 \cdot 10^{-3}$\\
			5.0 & $3.93 \cdot 10^{-3}$ & $ 3.93 \cdot 10^{-3}$  & $3.59 \cdot 10^{-3}$\\
			10.0 & $5.50 \cdot 10^{-4}$ & $5.50 \cdot 10^{-4}$ & $5.23 \cdot 10^{-4}$\\
			15.0 & $5.65 \cdot 10^{-5}$ & $5.65 \cdot 10^{-5}$ & $5.59 \cdot 10^{-5}$\\
			20.0 & $5.26 \cdot 10^{-6}$ & $5.26 \cdot 10^{-6}$ & $5.36 \cdot 10^{-6}$\\
			25.0 & $4.71 \cdot 10^{-7}$ & $4.71 \cdot 10^{-7}$ & $4.92 \cdot 10^{-7}$\\
			\hline
		\end{tabular}
	\end{center}
\end{table}



\begin{table}[!htb]
	\caption{ Calculated mean upcrossing rates for 50 eigenvalues
		\label{tab:6}}
	\begin{center}
		\begin{tabular}{c c c c}
			\hline
			$\eta = \zeta/ \lambda_1$   &  $\hat{\nu}_{Z_2}^+(\zeta) $  & $\check{\nu}_{Z_2}^+(\zeta) $ & $\tilde{\nu}_{Z_2}^+(\zeta)$ \\
			\hline
			2.0 & $6.55 \cdot 10^{-3}$ & $ 6.55 \cdot 10^{-3}$ & $ 5.93 \cdot 10^{-3}$\\
			5.0 & $3.25 \cdot 10^{-3}$ & $ 3.25 \cdot 10^{-3}$  & $2.98 \cdot 10^{-3}$\\
			10.0 & $5.04 \cdot 10^{-4}$ & $5.04 \cdot 10^{-4}$ & $4.70 \cdot 10^{-4}$\\
			15.0 & $5.44 \cdot 10^{-5}$ & $5.44 \cdot 10^{-5}$ & $5.28 \cdot 10^{-5}$\\
			20.0 & $5.19 \cdot 10^{-6}$ & $5.19 \cdot 10^{-6}$ & $5.20 \cdot 10^{-6}$\\
			25.0 & $4.71 \cdot 10^{-7}$ & $4.71 \cdot 10^{-7}$ & $4.86 \cdot 10^{-7}$\\
			\hline
		\end{tabular}
	\end{center}
\end{table}



To get a more detailed picture of how the crossing rate varies
with the number of eigenvalues retained, the mean upcrossing rate
was calculated for the level $\eta =20$ as a function of the
number of eigenvalues. The result is shown in
Figure~\ref{fig:crossingrates}. It was also decided to investigate
the effect of updating the truncated response representation so
that it had the correct variance. This was achieved by multiplying
the retained eigenvalues by a suitable factor. The effect of this
updating on the calculated upcrossing rate is also shown in
Figure~\ref{fig:crossingrates}. The figure indicates a couple of
interesting conclusions. Updating for variance can lead to
inaccurate results for the crossing rate for small to moderate
number of eigenvalues retained. Comparing
Figures~\ref{fig:eigenvalues} and \ref{fig:crossingrates} it is
seen that surprisingly accurate results are obtained for even a
small number of retained eigenvalues when the truncation is done
exactly where negative eigenvalues are followed by positive
eigenvalues. This seems to provide the right balance between the
terms in the response representation, and it indicates a useful
criterion for truncating the response representation for crossing
rate calculations.

\begin{figure}[!htb]
	\begin{center}
		\setlength{\unitlength}{0.012500in}%
		\dofig{6cm}{Figures/crossingrates.pdf}
	\end{center}
	\caption{The mean upcrossing rate of the level $\eta =20$ as a
		function of the number of retained eigenvalues.}
	\label{fig:crossingrates}
\end{figure}


It is also of great interest to observe that the simple Laplace
approximation in fact provides quite accurate estimates of the
mean upcrossing rates, and for this method the number of
eigenvalues has practically no effect on the computational burden.
Hence, from a practical point of view, this is an extremely
appealing method. In Table~\ref{tab:7}, the results
obtained by the hybrid method , the Laplace approximation and also
the simple approximation of Eq.~(\ref{nc15}) for 100 eigenvalues have been listed.
It is seen that while there is excellent agreement between the
hybrid method and the Laplace approximation, the simple
approximation leads to significantly lower values. In terms of
extreme value predictions, for the example structure at hand the
Laplace approximation is within about 1\% of the hybrid method,
while the simple approximations would lead to an underestimation
of typically 5-10\% compared with the two more accurate methods.


\begin{table}[!htb]
	\caption{ Calculated mean upcrossing rates for 100 eigenvalues
		\label{tab:7}}
	\begin{center}
		\begin{tabular}{c c c c}
			\hline
			$\eta = \zeta/ \lambda_1$   &  $\check{\nu}_{Z_2}^+(\zeta) $ & $\tilde{\nu}_{Z_2}^+(\zeta)$ & $\overline{\nu}_{Z_2}^+(\zeta)$ \\
			\hline
			2.0  & $6.17 \cdot 10^{-3}$ & $5.59 \cdot 10^{-3}$ & $6.03 \cdot 10^{-3}$\\
			5.0  & $3.03 \cdot 10^{-3}$ & $2.78 \cdot 10^{-3}$ & $2.74 \cdot 10^{-3}$\\
			10.0 & $4.71 \cdot 10^{-4}$ & $4.40 \cdot 10^{-4}$ & $3.65 \cdot 10^{-4}$\\
			15.0 & $5.11 \cdot 10^{-5}$ & $4.96 \cdot 10^{-5}$ & $3.44 \cdot 10^{-5}$\\
			20.0 & $4.88 \cdot 10^{-6}$ & $4.90 \cdot 10^{-6}$ & $2.94 \cdot 10^{-6}$\\
			25.0 & $4.44 \cdot 10^{-7}$ & $4.63 \cdot 10^{-7}$ & $2.44 \cdot 10^{-7}$\\
			\hline
		\end{tabular}
	\end{center}
\end{table}



\subsection{Wind excited structure}

The third example is a simple model of a wind excited structure. It could easily be adapted to the study of the response to wind loading of a high slender mast supporting a large light panel for illumination at a sports stadium. An example of such a structure is shown in Fig.~\ref{fig:lm}.


\begin{figure}[!htb]
	\begin{center}
		\setlength{\unitlength}{0.012500in}%
		\dofig{5cm}{Figures/stadium_mast-g.pdf}
	\end{center}
	\caption{A light-mast at a stadium.}
	\label{fig:lm}
\end{figure}

The theory can be extended to cover the general case of a MDOF
system as discussed by \cite{bepasp98}. However, for the purpose
of illustration, a SDOF model of structural response has been chosen.

Let the wind load on the structure be given as $F(t) = c\,
U(t)^2$, where $c$ is some constant, and $U(t)$ denotes a
stationary Gaussian process representing the wind speed. Writing
$U(t) = \bar{U}+V(t)$, where $\bar{U}$ = the mean wind speed, and
$V(t)$ denotes a zero mean Gaussian process, the
structural response to the loading process $\tilde{F}(t)= 2\,c\,
\bar{U}\, V(t) + c\,V(t)^2$ will be considered. That is, the constant force term $c\,
\bar{U}^2$ is neglected. For a linear dynamic model, its effect can be added at the end.

The structural response to $\tilde{F}(t)$ is assumed to be determined
by the equation,
\begin{equation} \label{wls:1}
   \ddot{Z}(t)+2\, \xi \,\omega_{0}\, \dot{Z}(t) +
   \omega_{0}^2 \,Z(t) = M^{-1} \tilde{F}(t)=\frac{2 c \bar{U}}{M} V(t) +\frac{c}{M} V(t)^2,
\end{equation}
where $\omega_{0}$ is the undamped natural frequency of the
system, $\xi$ is the relative damping and $M$ is the total mass.

The linear transfer function of this system is clearly given by
$\hat{L}(\omega)= [-\omega^2 \, +2\,i\,\xi\,\omega_0\,\omega\,+
\omega_{0}^2]^{-1}$. For this particular force model, the Volterra
series is of a degenerate kind. It can be shown that the QTF for
the quadratic response component of Eq.~\eqref{wls:1} is
given by the expression $c \, \hat{L}(\omega_i-\omega_j)/M$
\citep{naes87jsv}.

What is now needed to fully specify the $Q$-matrix of Eq.~\eqref{eq:cf12}
is the wind velocity spectrum $S_V(\omega)$. Here, a Davenport spectrum is adopted,
\begin{equation} \label{wa:1}
      \frac{\omega \,S_V(\omega)}{ \kappa \,\bar{U}^2} = \frac{4\, \theta^2}{(1
      +\,\theta^2)^{4/3}},
\end{equation}
where $\theta= \omega \, L/(2 \pi \bar{U})$, $\kappa$ is the
roughness/drag coefficient and $L$ is a length scale. For the
present example, a turbulence level of $\sigma /
\bar{U}=0.1$ and $0.3$ is assumed. The corresponding drag coefficient
$\kappa$ is given by the relation
$\sigma/\bar{U}=\sqrt{6\,\kappa}$. For the other parameters, $\bar{U}=50$ [m/s] and $L=1200$ [m] are selected. The damping ratio $\xi=0.1$ and the natural frequency $\omega_0=\pi$~[rad/sec].
Also, $c/M=0.5 \cdot 10^{-3}$. The frequency discretization of the
problem ranges from $-5$ to $5$~[rad/sec] with an increment of
$0.1$ [rad/sec].



In Table~\ref{tab:7.2} are listed calculated values for the mean upcrossing
rates $\nu_Z^+(\zeta)$ of the response process $Z(t)$ with no
restrictions imposed on the covariance structure for the case
$\sigma/\bar{U}=0.3$. To get an idea about the influence of the
correlation between $W$ and $\dot{W}$, the
upcrossing rate under the assumption that $\Sigma_{12}=0$ has also been calculated. This
crossing rate is denoted by $\tilde{\nu}_{Z}^{+}(\zeta)$. From
Table~\ref{tab:7.2}, it is seen that $\nu_Z^{+}(\zeta)$ and
$\tilde{\nu}_{Z}^{+}(\zeta)$ are almost equal, supporting the often
assumed negligible influence on the crossing rate of the
dependence between $W$ and $\dot{W}$.

Both estimates $\tilde{\nu}_Z^{+}(\zeta)$ and $\nu_{Z}^{+}(\zeta)$, are
then compared with the upcrossing rate $\nu_{g}^{+}(\zeta)$ of the
linear part of the response, which is a Gaussian process. The
spectrum of the linear part of the response is,
\begin{equation} \label{sz}
   S_Z(\omega) ={|\hat{L}(\omega)|}^2\,{\Big(\frac{2\,c\,\bar{U}}{M}\Big)}^2\,S_V(\omega).
\end{equation}


Comparing $\nu_{g}^{+}(\zeta)$ with $\nu_{Z}^{+}(\zeta)$, it is observed
that the relative difference between them increases as the level
$\zeta$ becomes higher, which is to be expected.



From a practical point of view, it is interesting to observe the
effect on the predicted extreme responses of the full quadratic
model for the wind loading as opposed to a Gaussian approximation.
For this purpose, the following extreme value distribution is
adopted,
\begin{equation}\label{}
F_{ext}(\zeta) = \exp \Big( - \nu_Z^+(\zeta)\,T \Big),
\end{equation}
where $F_{ext}(\zeta)$ is the cumulative distribution function of the
extreme response during a time interval of length $T$, that is, $max\{Z(t); 0 \le t \le T\}$.



Taking into account also the contributing constant forcing term
$c\,\bar{U}^2$ that was left out, and which adds a displacement of
0.13 m, it can be verified that the extreme values at the
exceedance probability level $10^{-2}$ are about 6-7\% higher for
the full quadratic model as compared to the Gaussian approximation
when $\sigma/\bar{U}=0.1$. The corresponding number is about 30\%
when $\sigma/\bar{U}=0.3$. Clearly, these numbers depend on the
specific examples considered, but they serve to indicate what the
effect will be of neglecting the quadratic part of the wind
loading.


\begin{table}[htb]
	\begin{center}
		\caption{Mean upcrossing rates against different levels.
			$\tilde{\nu}_{Z}^{+}(\zeta)$: numerical method assuming $\Sigma_{12}=0$;
			$\nu_{Z}^{+}(\zeta)$: numerical method; $\nu_{g}^{+}(\zeta)$:
			Gaussian approximation}\label{tab:7.2}
		\vspace{2mm}
		\begin{tabular}{c c c c c} 
			\hline
			{$\zeta~[\mbox{meters}]$} & {$\tilde{\nu}_Z^{+}(\zeta)$} &
			{$\nu_{Z}^{+}(\zeta)$} & {$\nu_{g}^{+}(\zeta)$} \\  
			\hline 
			0.1    & 0.29891  & 0.28714  &  0.22118  \\ 
			0.15   &  0.17232 & 0.16506  & 0.12021  \\ 
			0.2    &  0.08484 & 0.08092  & 0.05119  \\ 
			0.3    &  0.01454 & 0.01368  & 0.00446  \\ 
			0.4    &  0.00181 & 0.00167  & 1.46$\cdot 10^{-4}$  \\ 
			0.5    &  1.83$\cdot 10^{-4}$& 1.66$\cdot 10^{-4}$  & 1.82$\cdot 10^{-6}$  \\ 
			\hline
		\end{tabular}
	\end{center}
\end{table}









\section{Appendix 1 - The Average Crossing Rate}\label{Ap1}\index{average crossing rate}

A heuristic proof of Eqs.~\eqref{eq:cf19} and \eqref{eq:cf20} will be given. For a
rigorous derivation, the reader is referred to \cite{nar200}. The
derivation will be based on Parseval's formula\index{Parseval's formula}, which relates pairs of Fourier transforms. First, equation (18) will be rewritten. For
this purpose the Heaviside function\index{Heaviside function} $H(x)$ is introduced. It is
defined as follows: $H(x)=0 $ for $x<0$, $H(0)=1/2$ and $H(x)=1$ for
$x>0$. For simplicity, also write
$p(x)=f_{Z\dot{Z}}(\zeta, x)$, where $\zeta$ will be fixed. Then,
\begin{equation}\label{eq:cf51}
\nu_Z^+(\zeta) = \int_{-\, \infty}^\infty x\,p(x)\,H(x) \, dx.
\end{equation}

The Fourier transform $\mathcal{F}$ and its inverse
$\hat{\mathcal{F}}$ are (formally) defined as follows,
\begin{equation}\label{eq:cf52}
\hat{f}(t)= \mathcal{F}[f(x)] = \int_{-\, \infty}^\infty
f(x)\,e^{\ii \,x\,t}\, dx,
\end{equation}
and
\begin{equation}\label{eq:cf53}
f(x)= \hat{\mathcal{F}}[\hat{f}(t)] = \frac{1}{2\, \pi}\int_{-\,
\infty}^\infty \hat{f}(t)\,e^{-\,\ii \,x\,t}\, dt.
\end{equation}

Let $G(x)=x\,p(x)$. Parseval's formula provides us with the
following equation ($^*$ denotes complex conjugation),
\begin{equation}\label{eq:cf54}
\int_{-\, \infty}^\infty G(x)\, H(x) \,dx = \frac{1}{ 2\,\pi}
\int_{-\, \infty}^\infty \hat{G}(t)\,\hat{H}(t)^*\, dt.
\end{equation}
The Fourier transform $\hat{H}(t)$ of the Heaviside function can
be shown to be given by the relation \citep{brac86},
\begin{equation}\label{eq:cf55}
\hat{H}(t) = \frac{\ii}{t} + \pi \, \delta (t),
\end{equation}
where $\delta(\cdot)$ denotes Dirac's delta function\index{Dirac's delta function}. The Fourier transform pair $H(x)$, $\hat{H}(t)$ can only be interpreted as
such in a formal manner. A certain caution should therefore be
exercised when they are used in calculations. That is, one should
understand their limitations and proper use, which requires
knowledge about how the relation given in Eq.~\eqref{eq:cf55} was
established \citep{brac86}.

If $\hat{p}(t)$ denotes the Fourier transform of $p(x)$, then
$\hat{G}(t)= -\,\ii \, d \hat{p}(t)/dt$. From the definition of the
characteristic function, it is obtained that,
\begin{equation}\label{eq:cf56}
\hat{p}(t) =  \frac{1}{ 2\,\pi} \int_{-\, \infty}^\infty M_{Z
\dot{Z}}(s,t)\, e^{-\, \ii \zeta \,s} \, ds,
\end{equation}
leading to
\begin{equation}\label{eq:cf57}
\frac{d \hat{p}(t)}{dt} =  \frac{1}{ 2\,\pi} \int_{-\,
\infty}^\infty \frac{\partial M_{Z \dot{Z}}(s,t)}{\partial t}\,
e^{-\, \ii \zeta \,s} \, ds.
\end{equation}
Hence it is found that
\begin{eqnarray}\label{eq:cf58}
\nu_Z^+(\zeta) &=& \frac{-\,\ii}{(2\,\pi)^2}\int_{-\, \infty}^\infty
\int_{-\, \infty}^\infty \frac{\partial M_{Z
\dot{Z}}(s,t)}{\partial t}\, e^{-\, \ii \zeta \,s} \,
\bigg(\,-\,\frac{\ii}{t} + \pi\, \delta(t)\, \bigg) \, dt \,ds
\nonumber\\ &=& -\,\frac{1}{(2\,\pi)^2}\int_{-\, \infty}^\infty
\int_{-\, \infty}^\infty \frac{1}{t}\, \frac{\partial M_{Z
\dot{Z}}(s,t)}{\partial t}\, e^{-\, \ii \zeta \,s} \, dt \, ds
\nonumber\\ &-&\,\frac{\ii}{4\,\pi}\int_{-\, \infty}^\infty
\bigg(\frac{\partial M_{Z \dot{Z}}(s,t)}{\partial t}\bigg)_{t=0}
\, e^{-\, \ii \zeta \,s} \,ds,
\end{eqnarray}
which is seen to agree with Eq.~\eqref{eq:cf19}. A scrutiny of the
derivation of the (formal) Fourier transform of the Heaviside
function $H(\cdot)$, reveals that the integral with respect to $t$
in the above expression should be understood as a principal value
integral as described after Eq.~\eqref{eq:cf19}.

In the case of a stationary process,
$\nu_Z^+(\zeta)=\nu_Z(\zeta)/2$, where $\nu_Z(\zeta)$ denotes the
average crossing rate of the level $\zeta$. The corresponding Rice
formula for the crossing rate is given as,
\begin{equation}\label{eq:cf59}
\nu_Z(\zeta) = \int_{-\, \infty}^\infty x\,p(x)\,\sign(x) \, dx,
\end{equation}
where $\sign(x)=-\,1$ for $x<0$, $\sign(0)=0$, and $\sign(x)= 1$
for $x>0$. Since $\hat{\mathcal{F}}[\sign(x)]= 2\,\ii / t$, it is
seen from the derivation above that indeed,
\begin{eqnarray}\label{eq:cf60}
\nu_Z^+(\zeta) &=& \frac{-\,\ii}{(2\,\pi)^2}\int_{-\, \infty}^\infty
\int_{-\, \infty}^\infty \frac{\partial M_{Z
\dot{Z}}(s,t)}{\partial t}\, e^{-\, \ii \zeta \,s} \,
\bigg(\,-\,\frac{\ii}{t}\, \bigg) \, dt \,ds \nonumber\\ &=&
-\,\frac{1}{(2\,\pi)^2}\int_{-\, \infty}^\infty \int_{-\,
\infty}^\infty \frac{1}{t}\, \frac{\partial M_{Z
\dot{Z}}(s,t)}{\partial t}\, e^{-\, \ii \zeta \,s} \, dt \, ds.
\end{eqnarray}




\section{Appendix 2 - The Characteristic Function}\label{Ap2}\index{characteristic function}

\vspace{1mm} The derivation of the formula for the characteristic
function $M_{Z\dot{Z}}$ will be based on an integral equality
which is cited, but not proved, in a less general form by
\cite{cram46}. The prime (') will be used to denote transposition
of a matrix or vector, and $\det(A)$ denotes the determinant of
a square matrix $A$.

\begin{theorem}\label{1}
Let $A$ be a complex, symmetrical and nonsingular $n \times n$
matrix, and assume that the real part of $A$ is positive definite.
Let  $d$  denote a complex $n \times 1$ vector. Then the following
integral equality holds,
\begin{equation}\label{ie1}
\int_{-\infty}^\infty \cdots \int_{- \infty}^\infty \exp \bigg( d'
x - \frac{1}{2} \,\, x' \, A \, x \bigg) \, d x_1 \ldots d x_n =
\frac{(2\, \pi)^{n/2}}{\sqrt{\det(A)}}\, \exp \bigg(
\frac{1}{2}\,\, d'\, A^{-\,1} \, d \bigg).
\end{equation}
\end{theorem}

{\it Proof}: Let $\mu=A^{-\,1}\, d$. It is then easily verified
that $2\, d'\,x - x'\, A \, x = - \,(x-\mu)' \, A \, (x-\mu) + d'
\, A^{-\,1}\, d$. Hence, it follows that the desired result is
obtained if the following equality is proved ($z=x-\mu$),
\begin{equation}\label{ie2} 
\int_{-\infty - \ii \Im(\mu_1)}^{\infty - \ii \Im(\mu_1)} \cdots
\int_{- \infty - \ii \Im(\mu_n) }^{\infty - \ii \Im(\mu_n)} \exp
\bigg( - \frac{1}{2} \,\, z' \, A \, z \bigg) \, d z = \frac{(2\,
\pi)^{n/2}}{\sqrt{\det(A)}},
\end{equation}
where $\mu = ( \mu_1, \ldots, \mu_n)'$, $\Im(\mu_j)$ denotes the
imaginary part of $\mu_j$, and $dz=dz_1 \ldots dz_n$.

By assumption, $A= P + \ii Q$ where $P$ and $Q$ are real symmetric
matrices with $P$ positive definite. Then there exists a real
unitary matrix $U$ such that $U'PU=D= \diag(\lambda_1, \ldots,
\lambda_n)$, $\lambda_j > 0$ for $j=1, \ldots, n$. $\diag(\cdot)$
denotes a diagonal matrix with the indicated arguments on the
diagonal. The variable shift $z \to y$ is made, where $y=Uz$.
Then, $z'\, A \,z= y'D y + \ii \, y' \tilde{Q} y$, where
$\tilde{Q}=U'Q U$ is a real symmetric matrix. A scaling of the
coordinate variables is introduced by the variable shift $y \to
u$, where $y=\sqrt{D}u$, and $\sqrt{D}= \diag(\sqrt{\lambda_1},
\ldots, \sqrt{\lambda_n})$. Then, $z'\, A \,z= u'u + \ii \, u'
\hat{Q} u$, where $\hat{Q}=\sqrt{D}^{-\,1} \tilde{Q} \sqrt{D}^{-\,
1}$ is again a real symmetric matrix. Hence, a real unitary matrix
$V$ exists such that $V'\, \hat{Q}\,V= \Omega = \diag(\omega_1,
\ldots, \omega_n)$. Let $E = I + \ii \Omega = \diag(1+\ii \omega_1,
\ldots, 1+\ii \omega_n)$. By the variable shift $u \to v$, where $v=V
\, u$, it follows that $z'\, A \,z= v'\, E \, v = \sum_{j=1}^n
(1+\ii \omega_j)\,v_j^2$, $v = ( v_1, \ldots, v_n)'$. Since $U$ and
$V$ are unitary matrices, it is observed that
$dz=dy=(\sqrt{\det(D)})^{-\,1}\, du = (\sqrt{\det(D)})^{-\,1}\, dv
$. Denoting the integral on the left hand side of Eq.~(\ref{ie2}) by $I$, it is obtained that,
\begin{equation}\label{ie3}
I = \frac{1}{\sqrt{\det(D)}} \int_{-\infty - \ii c_1}^{\infty - \ii
c_1} \cdots \int_{- \infty - \ii c_n }^{\infty - \ii c_n} \exp \bigg(
- \frac{1}{2} \,\, \sum_{j=1}^n (1+\ii \omega_j)\,v_j^2 \bigg) \,
dv_1 \ldots dv_n,
\end{equation}
where $c_j$, $j=1, \ldots, n$ are suitable constants. Since the
functions $f_j(z)= \exp\{ - (1/2)\, (1+\ii \omega) \, z^2 \}$ are
entire functions in the complex variable $z=x+\ii y$, and since
$\lim_{|x| \to \infty} f(z) = 0$ uniformly in $y$ for $|y| \le
\mbox{const.}$, it follows that the path of integration from $-\,
\infty - \ii c_j$ to $\infty - \ii c_j$ in Eq.~(\ref{ie3}) can be
replaced by the path from $- \, \infty$ to $\infty$ by Cauchy's
theorem. From the standard result that $\int_{-\, \infty}^{\infty}
\exp \{ - a\, x^2 \} \, dx = \sqrt{\pi / a}$ for a complex
constant $a$ provided that $\Re(a)> 0$, it now follows that,
\begin{equation}\label{ie4}
I= \frac{1}{\sqrt{\det(D)}} \prod_{j=1}^n \frac{\sqrt{2\, \pi}}{
\sqrt{1+\ii \omega_j}}=\frac{1}{\sqrt{\det(D)}} \frac{(2\,
\pi)^{n/2}}{ \sqrt{\det(E)}}.
\end{equation}

To prove Eq.~(\ref{ie2}), it only remains to show that
$\det(D) \cdot \det(E) = \det(A)$. Invoking the fact that $U$ and
$V$ are unitary, it is obtained that,
\begin{eqnarray}\label{ie5}
\det(E)&=& \det(I + \ii \Omega)= \det(I+\ii \hat{Q}) = \det(I + \ii
\sqrt{D}^{-\, 1} \tilde{Q} \, \sqrt{D}^{-\, 1}) \nonumber\\ &=&
\det(\sqrt{D}^{-\, 1}\, D \, \sqrt{D}^{-\, 1} + \ii \sqrt{D}^{-\, 1}
\tilde{Q} \, \sqrt{D}^{-\, 1}) \nonumber\\ &=& (\det(D))^{-\, 1}
\, \det(D+\ii \tilde{Q})= (\det(D))^{-\, 1} \, \det(P+\ii Q) \nonumber\\
&=& \frac{\det(A)}{\det(D)},
\end{eqnarray}
which is what was needed to prove.


Our goal is to calculate the characteristic function,

\begin{equation} \label{charac:3}
M_{Z \dot{Z}}(u, v)=E \,[\, \exp \big( \ii u\, Z + \ii v\, \dot{Z}
\big) ].
\end{equation}


Using conditional probabilities, one may write,
\begin{equation} \label{charac:4}
M_{Z \dot{Z}}(u, v)= E \,\big[\, e^{\ii u\, Z} \,E \big[ e^{\ii v\,
\dot{Z}}|W \big] \big]= E \,\big[\, e^{\ii u\, Z} \,E \big[ e^{\ii v\,
(\dot{Z}|W)} \big] \big].
\end{equation}
From Eq.~\eqref{eq:cf17} it follows that $\dot{Z}= Y'\, \dot{W}$, where
$Y = \beta + 2 D W$, $\beta = (\beta_1, \ldots, \beta_n)'$ and $D=
\diag(\lambda_1, \ldots, \lambda_n)$. Hence, $(\dot{Z}|W)= Y' \,
(\dot{W}|W)$. $(\dot{W}|W=w)$, where $w=(w_1, \ldots, w_n)' \in
\mathbf{R}^n$, is a Gaussian vector with mean value $\mu=
\Sigma_{21}\,w $ and covariance matrix $C = \Sigma_{22} -
\Sigma_{21}\, \Sigma_{12}$, cf. \cite{ande58}. This means that
$(\dot{Z}|W=w)$ is a scalar Gaussian variable with mean value
$m=y'\, \mu= y'\, \Sigma_{21} \, w$ and variance $s^2= y'\, C \,
y$, where $y=\beta + 2 D w$. Invoking the expression for the
characteristic function of a Gaussian variable, it follows that
$E[\exp \{ \ii v\,(\dot{Z}|W)\}]= \exp \{ \ii v\, Y'\, \Sigma_{21}\,W
- \frac{1}{2}\, v^2 \, Y'\, C \, Y \}$. It is now obtained that,
\begin{equation} \label{charac:5}
M_{Z \dot{Z}}( u, v)= E \,[\, \exp \big( \ii u\, W' D W + \ii u \,
\beta' W + \ii v\, Y'\, \Sigma_{21}\,W - \frac{1}{2}\, v^2 \, Y'\, C
\, Y \big) ].
\end{equation}
The calculation of the expected value of this equation, amounts to
the calculation of the following integral,
\begin{equation}\label{}
M_{Z \dot{Z}}( u, v)= \frac{1}{(2 \pi)^{n/2}}\int_{-\infty}^\infty
\cdots \int_{- \infty}^\infty \exp \bigg( - \frac{1}{2}\, v^2\,
\beta'\,C \, \beta + d' w - \frac{1}{2} \,\, w' \, B \, w \bigg) \,
d w,
\end{equation}
where
\begin{equation}\label{}
d = (\ii \,u\, I + \ii \,v\,\Sigma_{12} - 2 \, v^2\, D \,C)\, \beta,
\end{equation}
and
\begin{equation}\label{}
B = I-2 \ii u\, D\, -4\,\ii \,v\,D \, \Sigma_{21}+ 4\,v^2 D \, C \,D.
\end{equation}

It is recognized by inspection that the matrix $B$ is not in
general symmetric. However, the expression $w'B w$ is a scalar
quantity, which implies that $w'B w = (w'B w)' = w'B'w$. Hence $w'B w = w' A w$, where $A$ is the symmetrized
version of $B$, that is, $A=(B+B')/2$. By invoking Theorem~\ref{1}, it
is therefore obtained:

\begin{theorem}\label{2}
Let the stochastic process $Z(t)$ be represented as given by
Eq.~\ref{eq:cf17}. Then the characteristic function $M_{Z \dot{Z}}(u,
v) = E \,[\, \exp \big( \ii u\, Z + \ii v\, \dot{Z} \big) ]$ of the
joint variable $(Z, \dot{Z})$ is given by the expression,
\begin{equation} \label{}
M_{Z \dot{Z}}( u, v)= \frac{1}{\sqrt{\det(A)}}\, \exp \bigg( -
\frac{1}{2}\, v^2\, \beta'\,C\, \beta + \frac{1}{2}\,\, d'\,
A^{-\,1} \, d \bigg),
\end{equation}
where
\begin{equation}\label{}
d = (\ii \,u\, I + \ii \,v\,\Sigma_{12} - 2 \, v^2\, D \,C)\, \beta,
\end{equation}
and
\begin{equation}\label{}
A = I-2\ii u\, D\, - 2\,\ii \,v\,\big( \,D \, \Sigma_{21}+
\Sigma_{12}\, D \big) + 4\,v^2 D \, C \,D.
\end{equation}
\end{theorem}

Note that in the absence of a linear component, that is, when
$Z(t) = \sum_{j=1}^n \lambda_j\, W_j(t)^2$, then,
\begin{equation} \label{}
M_{Z \dot{Z}}( u, v)= \frac{1}{\sqrt{\det(A)}}.
\end{equation}

The marginal characteristic functions $M_Z(u)= E[\exp\{\ii \, u\, Z
\}]$ and $M_{\dot{Z}}(v)=E[\exp\{\ii \, v\, \dot{Z} \}]$ are now
easily obtained by the relations $M_Z(u)=M_{Z \dot{Z}}( u, 0)$ and
$M_{\dot{Z}}(v)=M_{Z \dot{Z}}(0, v)$. In particular,
\begin{equation} \label{}
M_{Z }(u)= \frac{1}{\sqrt{\prod_{j=1}^n(1-2\,\ii \,u\, \lambda_j)}}\,
\exp \bigg( - \frac{1}{2}\sum_{j=1}^n \frac{ u^2\,
\beta_j^2}{1-2\, \ii \, u\, \lambda_j} \bigg),
\end{equation}
which is in agreement with previously derived results
\citep{kasi47,naes86aor}. It is also obtained that,
\begin{equation} \label{}
M_{\dot{Z}}(v)= \frac{1}{\sqrt{\det(\tilde{A})}}\, \exp \bigg( -
\frac{1}{2}\, v^2\, \beta'\,C\, \beta + \frac{1}{2}\,\, \tilde{d}'\,
\tilde{A}^{-\,1} \, \tilde{d} \bigg),
\end{equation}
where
\begin{equation}\label{}
\tilde{d} = ( \ii \,v\,\Sigma_{12} - 2 \, v^2\, D \,C)\, \beta,
\end{equation}
and
\begin{equation}\label{}
\tilde{A} = I - 2\,\ii \,v\,\big( \,D \, \Sigma_{21}+ \Sigma_{12}\, D
\big) + 4\,v^2 D \, C \,D.
\end{equation}



%\begin{figure}[!htb]
%\begin{center}
%\dofig{10.0cm}{Figures/windbsou01.pdf} \caption{The
%exceedance probability as a function of response level for
%$\sigma/\bar{U}=0.1$. $*\,*\,*$ : Full quadratic model;
%$\circ\,\circ\,\circ$ : Gaussian approximation.}
%\end{center}
%\end{figure}

%\begin{figure}[!htb]
%\begin{center}
%\dofig{10.0cm}{Figures/windasou03.pdf} \caption{The
%exceedance probability as a function of response level for
%$\sigma/\bar{U}=0.3$. $*\,*\,*$ : Full quadratic model;
%$\circ\,\circ\,\circ$ : Gaussian approximation.}
%\end{center}
%\end{figure}




