iccsa-21-wind

main.tex (34717B)

---

 \documentclass[runningheads]{llncs}
 
 \usepackage{amsmath}
 \usepackage{amssymb}
 \usepackage{booktabs}
 \usepackage{cite}
 \usepackage{graphicx}
 \usepackage{listings}
 \usepackage{url}
 \usepackage{textcomp}
 \usepackage{xcolor}
 
 \DeclareMathOperator{\Mode}{Mo}
 \DeclareMathOperator{\Expectation}{E}
 
 \begin{document}
 
 \title{Wind simulation using high-frequency velocity component measurements
     \thanks{Supported by Council for grants of the President of the Russian Federation (grant no.~\mbox{MK-383.2020.9}).}}
 \author{%
     Anton Gavrikov\orcidID{0000-0003-2128-8368} \and\\
     Ivan Gankevich\textsuperscript{*}\orcidID{0000-0001-7067-6928}
 }
 
 \lstdefinelanguage{rrr}{
     basicstyle={\rmfamily},
 	morekeywords={function},
     keywordstyle={\bfseries},
 	showstringspaces=false,
 	texcl=true,
 	escapeinside={LATEX}{END},
 	otherkeywords={function}
 }
 \lstset{
     language=rrr,
     literate={<-}{{$\gets$}}1,
 }
 
 \titlerunning{Wind simulation using high-frequency velocity component measurements}
 \authorrunning{A.\,Gavrikov, I.\,Gankevich}
 
 \institute{Saint Petersburg State University\\
     7-9 Universitetskaya Emb., St Petersburg 199034, Russia\\
     \email{st047437@student.spbu.ru},\\
     \email{i.gankevich@spbu.ru},\\
     \url{https://spbu.ru/}}
 
 \maketitle
 
 \begin{abstract}
 
     Wind simulation in the context of ship motions is used to estimate the
     effect of the wind on large containerships, sailboats and yachts. Wind
     models are typically based on a sum of harmonics with random phases and
     different amplitudes.  In this paper we propose to use autoregressive model
     to simulate the wind. This model is based on autocovariance function that
     can be estimated from the real-world data collected by anemometers.  We
     have found none of the data that meets our resolution requirements, and
     decided to produce the dataset ourselves using three-axis anemometer.  We
     built our own anemometer based on load cells, collected the data with the
     required resolution, verified the data using well-established statistical
     distributions, estimated autocovariance functions from the data and
     simulated the wind using autoregressive model. We have found that the load
     cell anemometer is capable of recording wind speed for statistical studies,
     but autoregressive model needs further calibration to reproduce the wind
     with the same statistical properties.
 
 \keywords{%
     load cell
     \and strain gauge
     \and anemometer
     \and three-dimensional ACF
     \and wind velocity PDF
     \and autoregressive model
     \and turbulence.
 }
 \end{abstract}
 
 \section{Introduction}
 
 Wind simulation in the context of ship motion simulation is the topic where
 multiple mathematical models are possible, and the choice of the model depends
 on the purpose of the model. In the course of the research where we apply
 autoregressive model to ocean wave simulation, we decided to investigate
 whether the same model can be used to simulate wind flow around ship hull.
 
 Wind simulation is studied at different scales and the closest scale for a ship
 in the ocean is wind turbine. One of the model that is used to
 describe air flow around wind turbines~\cite{veers1984modeling} is similar to
 Longuet---Higgins model which is typically used for ocean wave simulations.
 This model uses wind velocity frequency spectrum to determine coefficients
 and use them to generate wind velocity vector components.
 
 Wind velocity distribution is described by Kaimal
 spectrum~\cite{kaimal1972spectral,veers1988wind}
 \begin{equation*}
 S(f) = \frac {c_1 u_{*}^2 z / U(z)} {1 + c_2 \left(f z / U(z)\right)^{5/3} },\quad
 u_{*}=\frac{k U(z)}{\ln{(z/z_0)}},\quad c_1=105, \quad c_2=33,
 \end{equation*}
 Here \(u_{*}\) is shear velocity,
 \(k=0.4\) is von Karman constant, \(z_0\) is surface roughness
 coefficient, \(f\) is frequency,
 \(z\) is height above ground, \(U(z)\) is mean speed at height \(z\).
 
 This spectrum is used to simulate each wind velocity vector component at a
 specified point in space. For each component the same spectrum is used, but the
 coefficients are different~\cite{veers1984modeling}.
 %\begin{align*}
 %& c_1 = 12.3 && c_2 = 192 && \text{for }x,\\
 %& c_1 = 4    && c_2 = 70  && \text{for }y,\\
 %& c_1 = 0.5  && c_2 = 8   && \text{for }z.
 %\end{align*}
 The spectrum describes wind velocity vector in the plane that is perpendicular to
 the mean wind direction vector and travels in the same direction
 with mean wind speed. Time series is generated as Fourier series, coefficients
 of which are determined from the spectrum, and phases are random
 variables~\cite{veers1984modeling}:
 \begin{align*}
 & V(t) = \overline{V} + \sum\limits_{j=1}^{n}
 	\left( A_j \sin\omega_jt + B_j \cos\omega_jt\right), \\
 & A_j = \sqrt{\frac{1}{2} S_j \Delta\omega} \sin\phi_j, \quad
  B_j = \sqrt{\frac{1}{2} S_j \Delta\omega} \cos\phi_j.
 \end{align*}
 Here \(S_j\) is spectrum value at frequency \(\omega_j\),
 \(\phi_j\) is random variable which is uniformly distributed in \([0,2\pi]\).
 The result is one-dimensional vector-valued time series, each element of
 which is velocity vector at a specified point in time and space.
 
 In order to simulate wind velocity vector at multiple points in space,
 the authors use the function of coherency~--- the amount of correlation between
 wind speed at two points in space. This function has a form of an exponent and
 depends on frequency~\cite{veers1988wind}:
 \begin{equation*}
 \text{Coh}_{jk}(f) = \exp\left( -\frac {C\Delta r_{jk} f} {U(z)} \right),
 \end{equation*}
 where \(\Delta{}r_{jk}\) is the distance between \(i\) and \(j\) points and
 \(C\) is coherency decrement. Time series for each wind velocity vector
 component are generated independently and after that their spectra are modified
 in accordance with coherence function (the formulae are not presented here).
 
 In order to simulate wind with autoregressive model it is easier to use
 autocovariance function instead of spectra and coherence function. We can obtain
 autocovariance function from spectra as inverse Fourier transform using
 Wiener---Khinchin theorem. The formula that we obtained this way using
 various computer algebra programmes is too complex, but can be approximated by
 a decaying exponent:
 \begin{equation*}
 \gamma(t) = \sigma^2 \exp\left( -0.1 \frac{c_2^{3/5}}{c_1} t\right),
 \end{equation*}
 where \(\sigma^2\) is process variance (area under the spectrum).
 
 Unfortunately, this autocovariance function is one-dimensional and there is no
 easy way of obtaining three-dimensional analogue from the spectra.  In order
 to solve this problem we looked into datasets of wind speed measurements
 available for the research. However, most of them contain either one or two
 wind velocity vector components (in a form of wind speed and direction), they are difficult
 to get and their resolution is very small (we found only one paper that
 deals with three components~\cite{yim2000}). For our purposes we need
 resolution of at least one sample per second to simulate gusts and some form of
 turbulence. To summarise, our requirements for the dataset is to provide all
 three wind velocity vector components and has resolution of at least one sample per
 second.
 
 We failed to find such datasets and continued our research into anemometers
 that can generate the required dataset. One of the anemometers that can record
 all three components is ultrasonic
 anemometer~\cite{cosgrove2007ultra,arens2020anemometer,yakunin20173d}.  As
 commercial anemometers are too expensive for this research, we built our own
 version from the generally available electric components. However, this version
 failed to capture any meaningful data, and incidentally we decided to build an
 anemometer from load cells and strain gauges (which were originally intended for different research work). This anemometer is
 straightforward to construct, the electrical components are inexpensive, and it
 is easy to protect them from bad weather. This anemometer is able to record
 all three wind velocity vector components multiple times per second.
 
 In this paper we describe how load cell anemometer is built, then we collect
 dataset of all three wind speed components with one-second resolution, verify
 this anemometer using commercial analogue, verify measurements from the dataset
 using well-established distributions for wind speed and direction, and
 estimate autocovariance functions for autoregressive model from our dataset.
 Finally, we present preliminary wind simulation results using autoregressive model.
 
 \section{Methods}
 
 \subsection{Three-axis wind velocity measurements with load cell anemometer}
 
 In order to generate wind velocity field using four-dimensional (one temporal
 and three spatial dimensions) autoregressive model, we need to use
 four-dimensional wind velocity autocovariance function. Using Wiener---Khinchin
 theorem it is easy to compute the function from the spectrum. Unfortunately,
 most of the existing wind velocity historical data contains only wind velocity
 magnitude and direction. We can use them to reconstruct \(x\) and \(y\)
 spectrum, but there is no way to get spectrum for \(z\) coordinate from this
 data. Also, resolution of historical data is too small for
 wind simulation for ship motions. To solve these problems, we decided
 to build our own three-axis anemometer that measures wind velocity for all
 three axes at one point in space multiple times per second.
 
 To measure wind velocity we used resistive foil strain gauges mounted on the
 arms aligned perpendicular to the axes directions. Inside each arm we placed
 aluminium load cell with two strain gauges: one on the bending side and another
 one on the lateral side. Load cells use Wheatstone half-bridge to measure the
 resistance of the gauges and are connected to the circuit that measures the
 resistance and transmits it to the microcontroller in digital format.
 Microcontroller then records the value for each load cell and transmits all of
 them in textual form to the main computer.  The main computer then adds a
 timestamp and saves it to the database.  3-D anemometer model is shown in
 figure~\ref{fig-anemometer}.
 
 \begin{figure}
     \centering
     \begin{minipage}{0.55\textwidth}
         \includegraphics{build/inkscape/anemometer.eps}
     \end{minipage}
     \begin{minipage}{0.35\textwidth}
         Anemometer properties:\phantom{12345}
         \begin{tabular}{ll}
             \toprule
             Load cell capacity & 1 kg \\
             Load cell amplifier & HX711 \\
             Microcontroller & ATmega328P \\
             \bottomrule
         \end{tabular}
     \end{minipage}
     %\includegraphics{build/gnuplot/anemometer.eps}
     \caption{Three-axis anemometer. Arms are inserted into the housing
     and fixed using bolts that go through the circular holes. Red, green, blue colours
 denote load cells for \(x\), \(y\), \(z\) axes respectively.\label{fig-anemometer}}
 \end{figure}
 
 Each load cell faces the direction that is perpendicular to the directions of
 other load cells. When the wind blows in the direction of particular load cell, only
 this cell bends. When the wind blows in an arbitrary direction which is not
 perpendicular to any load cell faces, then all load cells bend, but the
 pressure force is smaller. Pressures of all load cells are recorded
 simultaneously, and we can use Bernoulli's equation to compute wind velocity
 from them.
 
 Bernoulli's equation is written as
 \begin{equation}
     \rho\frac{\upsilon^2}{2} + \rho gz = p_0-p.
 \end{equation}
 Here \(\upsilon\) is wind velocity magnitude, \(g\) is gravitational acceleration,
 \(z\) is vertical coordinate, \(\rho\) is air density, \(p\) is pressure on the
 load cell and \(p_0\) is atmospheric pressure. Pressure force \(\vec{F}\) acting on a load
 cell is written as
 \begin{equation}
     \vec{F} = p S \vec{n},
 \end{equation}
 where \(S\) is area of the side of the load cell on which the force is applied
 and \(\vec{n}\) is normal vector. Arms in which load cells reside have 
 front and back side with much larger areas than the left and right side, therefore
 we neglect forces on them. Additionally, on the ground \(\rho{}gz\) term vanishes.
 
 For all load cells we have a system of three equations
 \begin{equation}
     \begin{cases}
         \upsilon_x^2 \propto F_x \\
         \upsilon_y^2 \propto F_y \\
         \upsilon_z^2 \propto F_z \\
     \end{cases}.
 \end{equation}
 Hence \(\upsilon_{x}=\alpha_{x}\sqrt{F_{x}}\),
 \(\upsilon_{y}=\alpha_{y}\sqrt{F_{y}}\),
 \(\upsilon_{z}=\alpha_{z}\sqrt{F_{z}}\) where \(\alpha_{x,y,z}\) are constants of
 proportionality. Therefore, to obtain wind velocity we take \emph{square root} of the 
 value measured by the load cell and multiply it by some coefficient determined
 empirically during anemometer calibration.
 
 \subsection{Per-axis probability distribution function for wind velocity}
 
 Scalar wind velocity is described by Weibull distribution~\cite{justus1976}.
 Weibull probability distribution function is written as
 \begin{equation}
     f\left(\upsilon;b,c\right) = b c \left(b \upsilon\right)^{c-1}
     \exp\left(-\left(b\upsilon\right)^{c}\right).
 \end{equation}
 Here \(\upsilon>0\) is scalar wind velocity, \(b>0\) is scale parameter and \(c>0\)
 is shape parameter. This function is defined for positive wind velocity, since
 scalar wind velocity is a length of wind velocity vector
 \(\upsilon=\sqrt{\smash[b]{\upsilon_x^2+\upsilon_y^2+\upsilon_z^2}}\). However, projection of
 wind velocity vector on \(x\), \(y\) or \(z\) axis may be negative.
 Our solution to this problem is to use two Weibull distributions: one for positive and
 one for negative projection~--- but with different parameters.
 \begin{equation}
     f\left(\upsilon_x;b_1,c_1,b_2,c_2\right) = \begin{cases}
         b_1 c_1 \left(b_1 \left|\upsilon_x\right|\right)^{c_1-1} \exp\left(-\left(b_1\left|\upsilon_x\right|\right)^{c_1}\right).
         \qquad \text{if } \upsilon_x < 0 \\
         b_2 c_2 \left(b_2 \left|\upsilon_x\right|\right)^{c_2-1} \exp\left(-\left(b_2\left|\upsilon_x\right|\right)^{c_2}\right).
         \qquad \text{if } \upsilon_x \geq 0 \\
 
     \end{cases}\label{eq-velocity-distribution}.
 \end{equation}
 Here \(b_{1,2}>0\) and \(c_{1,2}>0\) are parameters of the distribution
 that control the scale and the shape, \(\upsilon_x\) is the projection of the velocity
 vector on \(x\) axis. The same formula is used for \(y\) and \(z\) axis.
 
 %We could simply
 %extend Weibull distribution for negative values, but this may lead to multimodal distribution
 %because the maximum value may not be located at \(\upsilon=0\). Our solution to
 %this problem is to make substitution \(\upsilon_x'=\left|\upsilon_x-\Mode\upsilon_x\right|\),
 %where \(\Mode\upsilon_x\) is the mode of wind velocity projection to the corresponding axis,
 %and modulus makes the graph symmetric with respect to \(\upsilon_x'=0\).
 %The substitution yields
 %\begin{equation}
 %    f\left(\upsilon_x';b,c\right) = b c
 %    \exp\left(-\left|b\upsilon_x'-g\left(c\right)\right|^{c}\right),
 %    \qquad
 %    g\left(c\right) = \left(\frac{c-1}{c}\right)^{1/c},
 %\end{equation}
 %which is the formula for the right side of the distribution. In order to get the
 %formula for the left side, we have to use an additional set of parameters \(b\), \(c\).
 %Then the final formula for both sides is written as
 %\begin{equation}
 %    f\left(\upsilon_x';b_1,c_1,b_2,c_2\right) = \begin{cases}
 %        b_1 c_1 \exp\left(-\left|b_1\upsilon_x'-g\left(c_1\right)\right|^{c_1}\right)
 %        \qquad \text{if } \upsilon_x < \Mode\upsilon_x \\
 %        b_2 c_2 \exp\left(-\left|b_2\upsilon_x'-g\left(c_2\right)\right|^{c_2}\right)
 %        \qquad \text{if } \upsilon_x \geq \Mode\upsilon_x \\
 %
 %    \end{cases}\label{eq-velocity-distribution}.
 %\end{equation}
 %Here \(b_{1,2}>0\) and \(c_{1,2}>0\) are parameters of the distribution
 %that control the scale and the shape.
 
 \subsection{Three-dimensional ACF of wind velocity}
 
 Usually, autocovariance is modelled using exponential functions~\cite{box1976time}.
 In this paper we use one-dimensional autocovariance function written as
 \begin{equation}
     K\left(t\right) = a_3 \exp\left(-\left(b_3 t\right)^{c_3} \right).
     \label{eq-acf-approximation}
 \end{equation}
 Here \(a_3>0\), \(b_3>0\) and \(c_3>0\) are parameters of the autocovariance
 function that control the shape of the exponent.
 
 In order to construct three-dimensional autocovariance function we assume that
 one-dimensional autocovariance function is the same for each coordinate and
 multiply them.
 \begin{equation}
     K\left(t,x,y,z\right) = a \exp\left(
         -\left(b_t t\right)^{c_t}
         -\left(b_x x\right)^{c_x}
         -\left(b_y y\right)^{c_y}
         -\left(b_z z\right)^{c_z}
     \right).
     \label{eq-acf}
 \end{equation}
 Here \(a>0\), \(b_{t,x,y,z}>0\) and \(c_{t,x,y,z}>0\) are parameters of the autocovariance
 function. Parameter \(a\) and \(\exp\left(b\right)\) are proportional to wind
 velocity projection on the corresponding axis.
 Parameter \(c\) controls the shape of the autocovariance function in
 the corresponding direction; it does not have simple relationship to the wind
 velocity statistical parameters.
 
 \subsection{Data collection and preprocessing}
 
 We installed anemometer on the tripod and placed it on the balcony. Then we
 connected load cells to the microcontroller via HX711 load cell amplifiers and
 programmed the microcontroller to record the output of each sensor every second
 and print it on the standard output. Then we connected the microcontroller to
 the computer via USB interface and wrote a script to collect the data coming
 from the USB and store it in the SQLite database. We decided to store raw
 sensor values in the range from 0 to 65535 to be able to calibrate anemometer
 later.
 
 We calibrated three-axis anemometer using commercial anemometer HP-866A and a
 fan that rotates with constant speed. First, we measured the wind speed that
 the fan generates when attached to commercial anemometer. Then we successively
 placed the fan behind and in front of each arm of our anemometer and measured
 values that the corresponding load cell reported for each side of the arm with
 and without the fan. Then we were able to calculate two coefficients for each
 axis: one for negative and another one for positive wind velocities. The
 coefficient equals the raw sensor value that is equivalent to the wind speed of 1~m/s
 (table~\ref{tab-coefficients}).
 
 \begin{table}
     \begin{minipage}[t]{0.35\textwidth}
         \centering
         \begin{tabular}{lll}
             \toprule
             Axis & \(C_1\) & \(C_2\) \\
             \midrule
             X & 11.19 & 12.31 \\
             Y & 11.46 & 11.25 \\
             Z & 13.55 & 13.90 \\
             \bottomrule
         \end{tabular}
         \caption{Calibration coefficients for each arm of three-axis anemometer:
             \(C_1\) is for negative values and \(C_2\) is for positive values.
         \label{tab-coefficients}}
     \end{minipage}
     ~
     \begin{minipage}[t]{0.55\textwidth}
         \centering
         \begin{tabular}{ll}
             \toprule
             Time span & 36 days \\
             Size & 122 Mb \\
             No. of samples & 3\,157\,234 \\
             No. of samples after filtering & 2\,775\,387 \\
             Resolution & 1 sample per second \\
             \bottomrule
         \end{tabular}
         \caption{Dataset properties.\label{tab-dataset}}
     \end{minipage}
 \end{table}
 
 We noticed that ambient temperature affects values reported by our load cells:
 when the load cell heats up (cools down), it reports values that increase
 (decrease) linearly in time due to thermal expansion of the material. We
 removed this linear trend from the measured values using linear regression. The
 code in R~\cite{r-language} that transforms raw sensor values into wind speed is
 presented in listing~\ref{lst-sample-to-speed}.
 
 \begin{lstlisting}[label={lst-sample-to-speed},caption={The code that transforms raw load cell sensor values into wind speed projections to the corresponding axis.}]
 sampleToSpeed <- function(x, c1, c2) {
  t <- c(1:length(x))
  reg <- lm(x~t)
  x <- x - reg$fitted.values LATEX{\color{black!70}\textit{\# remove linear trend}}END
  x <- sign(x)*sqrt(abs(x))  LATEX{\color{black!70}\textit{\# convert from force to velocity}}END
  x[x<0] = x[x<0] / c1      LATEX{\color{black!70}\textit{\# scale sensor values to wind speed}}END
  x[x>0] = x[x>0] / c2      LATEX{\color{black!70}\textit{\# using calibration coefficients}}END
  x
 }
 \end{lstlisting}
 
 Over a period of one month we collected 3.1M samples and filtered out 12\% of
 them because they had too large unnatural values. We attributed these values to
 measurement errors as they spread uniformly across all the time span and are
 surrounded by the values of regular magnitude.
 %From the remaining data we choose days with wind speeds above the average as
 %reported by EMERCOM of Russia\footnote{https://en.mchs.gov.ru/} for Saint
 %Petersburg.
 After that we divided each day into two-hour intervals over which
 we collected the statistics individually.
 %From these intervals we choose the
 %ones with data distributions close to unimodal, so we could easily fit them
 %into the model.
 The statistics for each interval is presented in
 figure~\ref{fig-intervals}, dataset properties are presented in table~\ref{tab-dataset}.
 
 Unique feature of three-axis anemometer is that it measures both velocity of
 incident air flow towards the arms and the turbulent flow that forms behind the
 arms. Turbulent flow velocity distribution peak is often smaller than the peak
 of incident flow. To exclude it from the measurements one can choose the side
 with the largest peak (either positive or negative part
 of~\eqref{eq-velocity-distribution}), but in this paper we left them as is
 for the purpose of the research.
 
 \begin{figure}
     \centering
     \includegraphics{build/daily-stats.eps}
     \caption{The statistics for each interval of the collected data. Here
         \(\upsilon_x\), \(\upsilon_y\) and \(\upsilon_z\) are mean speeds for
         each interval for the corresponding axes, \(\upsilon\) is mean scalar wind speed
         for each interval. The horizontal line shows overall average speed.
         Yellow rectangles denote days when EMERCOM of Russia (https://en.mchs.gov.ru/)
         reported wind speeds above average.
         \label{fig-intervals}}
 \end{figure}
 
 \section{Results}
 
 \subsection{Anemometer verification}
 
 In order to verify that our anemometer produces correct measurements we
 calculated wind speed and direction from the collected samples and fitted them
 into Weibull distribution and von Mises distribution respectively.  These are
 typical models for wind speed and direction~\cite{carta2008,carta2008joint}.
 Then we found the intervals with the best and the worst fit for these models
 using normalised root-mean-square error (NRMSE) calculated as
 \begin{equation}
     \text{NRMSE} = \frac{\sqrt{\Expectation\left[\left(X_\text{observed}-X_\text{estimated}\right)^2\right]}}{X_\text{max}-X_\text{min}}.
 \end{equation}
 Here \(\Expectation\) is statistical mean, \(X_\text{observed}\) and
 \(X_\text{estimated}\) are observed and estimated values respectively.
 
 The wind speed data collected with three-axis anemometer was approximated by
 Weibull distribution using least-squares fitting. Negative and positive wind
 speed projections to each axis both have this distribution, but with different
 parameters. Most of the data intervals contain only one prevalent mean wind
 direction, which means that one of the distributions is for incident wind flow
 on the arm of the anemometer and another one is for the turbulent flow that
 forms behind the arm.  For \(z\) axis both left and right distributions have
 similar shapes, for \(x\) and \(y\) axes the distribution for incident flow is
 taller than the distribution for turbulent flow. The best-fit and worst-fit
 distributions for each axis are presented in figure~\ref{fig-velocity-distributions}.
 
 \begin{figure}
     \centering
     \includegraphics{build/gnuplot/velocity-dist.eps}
     \caption{Per-axis wind velocity distributions fitted into Weibull
         distribution.
         Data intervals with the largest error (the first row),
         data intervals with the smallest error (the second row). Red line shows estimated
         probability density of positive wind speed projections, blue line shows estimated
         probability density of negative wind speed projections and circles denote observed
         probability density of wind speed projections.
     \label{fig-velocity-distributions}}
 \end{figure}
 
 Wind direction was approximated by von Mises distribution using least-squares
 fitting. Following the common
 practice~\cite{feng2015sectors} we divided direction axis
 into sectors: from -180\textdegree{} to -90\textdegree{}, from -90\textdegree{} to
 0\textdegree{}, from 0 to 90\textdegree{} and from 90\textdegree{} to
 180\textdegree{}~--- and fitted each sector independently. We chose four sector
 to have one sector for each side of the anemometer.
 The best-fit and worst-fit distributions for each sector
 are presented in figure~\ref{fig-direction-distributions}.
 
 \begin{figure}
     \centering
     \begin{minipage}{0.49\textwidth}
         \centering
         \includegraphics{build/gnuplot/velocity-xy-dist-max.eps}
     \end{minipage}
     \begin{minipage}{0.49\textwidth}
         \centering
         \includegraphics{build/gnuplot/velocity-xy-dist-min.eps}
     \end{minipage}
     \vskip2\baselineskip
     \includegraphics{build/gnuplot/direction-dist.eps}
     \caption{Bivariate \(x\) and \(y\) velocity distribution
         in polar coordinates (first row), wind direction distributions fitted into
         von Mises distribution (second row). Data interval with the largest error (left column), data
         interval with the smallest error (right column).
         Red line shows estimated
         probability density of positive direction angles, blue line shows estimated
         probability density of negative direction angles and circles denote observed
         probability density of direction angles. 0\textdegree{} is north.
     \label{fig-direction-distributions}}
 \end{figure}
 
 In order to verify that our anemometer produces correct time series we
 performed synchronous measurements with commercial anemometer HP-866A. Two
 anemometers were placed in the open field near the shore. The commercial
 anemometer was directed towards the mean wind direction as it measures the speed
 only. The data from the two anemometers was recorded synchronously. From the
 data we selected intervals with wind gusts and compared the measurements of the
 two anemometers. To compare them we scaled load cell anemometer measurements to
 minimise NRMSE of the difference between the two time series to compensate for
 the errors in calibration coefficients. The results showed that there is some
 correspondence between the measurements of the two anemometers
 (figure~\ref{fig-hold-peak}).
 
 \begin{figure}
     \centering
     \includegraphics{build/gnuplot/hold-peak.eps}
     \caption{Comparison of measurements obtained from three-axis and HP-866A
     anemometers.\label{fig-hold-peak}}
 \end{figure}
 
 Finally, we computed autocovariance for each axis as
 \begin{equation}
     K(\tau) = \Expectation\left[ \left(X_t-\bar{X}\right)\left(X_{t-\tau}-\bar{X}\right) \right]
 \end{equation}
 and fitted it into~\eqref{eq-acf-approximation}.
 Per-axis ACFs have pronounced peak at nought lag and long tails.
 The largest NRMSE is 2.3\%.
 Variances for \(x\) and \(y\) axes are comparable, but ACF for \(z\) axis has
 much lower variance.
 Parameters \(a\) and \(b\) from~\eqref{eq-acf-approximation} are positively
 correlated with wind speed for the corresponding axis.
 The best-fit and worst-fit ACFs for each axis
 are presented in figure~\ref{fig-acf}.
 
 \begin{figure}
     \centering
     \includegraphics{build/gnuplot/acf.eps}
     \caption{Per-axis wind velocity ACF fitted into~\eqref{eq-acf-approximation}.
         Data intervals with the largest error (the first row),
         data intervals with the smallest error (the second row). Red line shows estimated
         ACF of positive wind speed projections, blue line shows estimated
         ACF of negative wind speed projections and circles denote observed
         ACF of wind speed projections. Note: graphs use mirrored axes.
     \label{fig-acf}}
 \end{figure}
 
 \subsection{Turbulence coefficient}
 
 Since our anemometer measures both incident and turbulent flow we took an
 opportunity to study the speed of the turbulent flow in relation to incident
 flow. We calculated the absolute mean values of positive and negative wind
 velocity projections for each time interval of the dataset. We considered the
 flow with the larger absolute mean value incident, and the other one is
 turbulent.  Then we calculated turbulence coefficient as the ratio of absolute
 mean speed of turbulent flow to the absolute mean speed of incident flow. We
 found that the average ratio is close to 60\%. It seems, that the ratio
 decreases as the wind speed increases for wind speeds of 1--5 m/s, but for
 higher wind speeds we do not have the data. Turbulence coefficient
 can be used in wind simulation models to control the magnitude of turbulent
 flow that forms behind the obstacle~\cite{gavrikov2020wind}.
 
 \begin{figure}
     \centering
     \includegraphics{build/gnuplot/turbulence.eps}
     \caption{The ratio of absolute mean speed of turbulent flow to the absolute mean
     speed of incident flow for each axis.\label{fig-turbulence}}
 \end{figure}
 
 \subsection{Wind simulation using measured ACFs}
 
 We simulated three-dimensional wind velocity using autoregressive model
 implemented in Virtual Testbed~\cite{bogdanov2020vtestbed}~--- a
 programme for workstations that simulates ship motions in extreme conditions
 and physical phenomena that causes them (ocean waves, wind, compartment
 flooding etc.).  Using the data obtained with three-axis anemometer on March 28,
 2021, 01:00-03:00 UTC we approximated four-dimensional autocovariance function
 using~\eqref{eq-acf} by setting the corresponding parameters from
 one-dimensional autocovariance functions estimated from the data obtained with
 anemometer, all other parameters were set close to nought. Non-nought
 parameters are listed in table~\ref{tab-arma-parameters}.  We found that
 velocity and direction distributions and ACFs of each axis of simulated wind
 and real wind are similar in shape, but are too far away from each other
 (see~figure~\ref{fig-vtestbed}).  We consider these results preliminary and
 will investigate them further in future work.
 
 \begin{table}
     \centering
     \caption{Input parameters for AR model that were used to simulate wind velocity.
     \label{tab-arma-parameters}}
     \begin{tabular}{lllll}
         \toprule
         Axis & ACF \(a\) & \phantom{0}ACF \(b\) & \phantom{0}ACF \(c\) & \phantom{0}Mean velocity, m/s \\
         \midrule
         \(x\) & 1.793 & \phantom{0}0.0214 & \phantom{0}0.2603 & \phantom{0}-2.439 \\
         \(y\) & 1.423 & \phantom{0}0.01429 & \phantom{0}0.2852 & \phantom{0}-2.158 \\
         \(z\) & 0.9075 & \phantom{0}0.06322 & \phantom{0}0.3349 & \phantom{0}-1.367 \\
         \bottomrule
     \end{tabular}
 \end{table}
 
 \begin{figure}
     \centering
     \includegraphics{build/gnuplot/vtestbed.eps}
     \caption{Comparison of simulated wind velocity and direction distributions
         and per-axis ACFs
     between the data from the anemometer and the data from Virtual Testbed
     (2021, March 28, 01:00--03:00, UTC).\label{fig-vtestbed}}
 \end{figure}
 
 \section{Discussion}
 
 NRMSE of wind speed distribution approximation has positive correlation with
 wind speed: the larger the wind speed, the larger the error and vice versa.
 Larger error for low wind speeds is caused by larger skewness and kurtosis (see
 the first row of figure~\ref{fig-velocity-distributions}). Similar
 approximation errors can be found in~\cite{carta2008joint} where the authors
 improve approximation accuracy using joint wind speed and direction
 distributions. Such studies are outside of the scope of this paper, because
 here we verify anemometer measurements using well-established mathematical
 models, but the future work may include the study of these improvements.
 
 NRMSE of wind direction distribution approximation has negative correlation with
 wind speed: the larger the wind speed, the smaller the error and vice versa.
 This is in agreement with physical laws: the faster the flow is the more
 determinate its mean direction becomes, and the slower the flow is the more
 indeterminate its mean direction is.
 
 Three-axis anemometer disadvantages are the following.
 The arm for the \(z\) axis is horizontal, and snow and rain put additional load
 on this cell distorting the measurements. 
 Also, thermal expansion and contraction of the material changes
 the resistance of load cells and distorts the measurements. 
 Pressure force on the arm is exerted by individual air particles and
     is represented by choppy time series, as opposed to real physical signal
     that is represented by smooth graph.
 The first two defficiences can be compensated in software by removing linear trend
 from the corresponding interval. The last one make anemometer useful only for
 offline studies, i.e.~it is useful to gather statistics, but is unable to measure
 immediate wind speed and direction.
 
 We used a balcony for long-term measurements and open field for verification
 and calibration. We found no clues that the balcony affected the distributions
 and ACFs of wind speed. The only visible effect is that the wind direction is
 always parallel to the wall which agrees with physical laws. Since we measure
 pressure force directly, the mean wind direction does not affect the form
 of the distributions, but only their parameters.
 
 \section{Conclusion}
 
 In this paper we proposed three-axis anemometer that measures wind speed for
 each axis independently. We analysed the data collected by this anemometer and
 verified that per-axis wind speeds fit into Weibull distribution with the
 largest NRMSE of 6.7\% and wind directions fit into von Mises distribution with
 the largest NRMSE of 11\%. We estimated autocovariance functions for wind speed
 for each axis of the anemometer and used this approximations to simulate wind
 flow in Virtual Testbed. The parameters of these functions allow to control
 both wind speed and mean direction. The future work is to construct an array of
 anemometers that is able to measure spatial autocovariance using the proposed
 anemometer as the base.
 
 \subsubsection*{Acknowledgements.}
 Research work is supported by Council for grants of the President of the
 Russian Federation (grant no.~MK-383.2020.9).
 
 \bibliographystyle{splncs04}
 \bibliography{main.bib}
 
 \end{document}