Consider historical observation dates $t_0 < t_1 < \cdots < t_n$.
From the state variable equations
\begin{align*}
dZ_t^1&=-kZ_t^1dt+h_1dW_t^1+h_2dW_t^2,\\
dZ_t^2&=h_0dW_t^1.
\end{align*}
We obtain that, for $i=1, \ldots, n$,
\begin{align*}
Z_{t_i}^1 &= e^{-k \Delta t_i} Z_{t_{i-1}}^1 + h_1 \int_{t_{i-1}}^{t_i}e^{-k (t_i-s)}dW_s^1 + h_2 \int_{t_{i-1}}^{t_i}e^{-k (t_i-s)}dW_s^2,\tag{1}\\
Z_{t_i}^2 &= Z_{t_{i-1}}^2 + h_0 \int_{t_{i-1}}^{t_i}dW_s^1.\tag{2}
\end{align*}
Let
\begin{align*}
\pmb{x}_{t_i} &= [Z_{t_i}^1, \ Z_{t_i}^2]^T,
\end{align*}
and
\begin{align*}
F &= \left(\!
\begin{array}{cc}
e^{-k\Delta t_i} & 0\\
0 & 1
\end{array}
\!\right).
\end{align*}
Moreover, let
\begin{align*}
\pmb{v}_{t_i} = \bigg[h_1 \int_{t_{i-1}}^{t_i}e^{-k (t_i-s)}dW_s^1 + h_2 \int_{t_{i-1}}^{t_i}e^{-k (t_i-s)}dW_s^2, \ h_0 \int_{t_{i-1}}^{t_i}dW_s^1\bigg]^T
\end{align*}
be a two dimensional normal random vector with zero mean and covariance matrix
\begin{align*}
Q = \left(\!
\begin{array}{cc}
\frac{h_1^2+h_2^2}{2k}\big(1-e^{-2k \Delta t_i} \big) & \frac{h_0h_1}{k} \big(1-e^{-k\Delta t_i} \big) \\
\frac{h_0h_1}{k} \big(1-e^{-k\Delta t_i} \big) & h_0^2 \Delta t_i
\end{array}
\!\right).
\end{align*}
Then, based on $(1)$ and $(2)$,
\begin{align*}
\pmb{x}_{t_i} = F \pmb{x}_{t_{i-1}} + \pmb{v}_{t_i},
\end{align*}
which is the transition equation or state equation.
Let $T_1 < \cdots < T_m$ be the futures maturities. From the equation
\begin{align*}
\ln F(t,T)&=\ln F(0,T)+\left(Z_t^1 e^{-k(T-t)}+Z_t^2\right) \\
&\qquad -\frac{1}{4k}\left[(1-e^{-2kT})(h_1^2+h_2^2))+4h_1h_0(1-e^{-kT})+2h_0^2tk\right],
\end{align*}
we obtain that
\begin{align*}
\ln F(t_i,T_j)&=\ln F(0,T_j)+\left(Z_{t_i}^1 e^{-k(T_j-t_i)}+Z_{t_i}^2\right) \\
&\qquad -\frac{1}{4k}\left[(1-e^{-2kT_j})(h_1^2+h_2^2))+4h_1h_0(1-e^{-kT_j})+2h_0^2t_ik\right].\tag{3}
\end{align*}
For observation time $t_i$, let $\pmb{y}_{t_i}$ be an $m$-dimensional observation vector defined by
\begin{align*}
\pmb{y}_{t_i} =
\left(\!
\begin{array}{c}
\ln F(t_i, T_1)\\
\vdots\\
\ln F(t_i, T_m)
\end{array}
\!\right),
\end{align*}
and $\pmb{d}_{t_i}$ be an $m$-dimensional deterministic vector defined by
\begin{align*}
\pmb{d}_{t_i} =
\left(\!
\begin{array}{c}
\ln F(0,T_1) - \frac{1}{4k}\left[(1-e^{-2kT_1})(h_1^2+h_2^2))+4h_1h_0(1-e^{-kT_1})+2h_0^2t_ik\right]\\
\phantom{\frac{\frac{1}{1}}{\frac{1}{1}}} \vdots \phantom{\frac{\frac{1}{1}}{\frac{1}{1}}}\\
\ln F(0,T_m) - \frac{1}{4k}\left[(1-e^{-2kT_m})(h_1^2+h_2^2))+4h_1h_0(1-e^{-kT_m})+2h_0^2t_ik\right]
\end{array}
\!\right).
\end{align*}
Moreover, let $H$ be an $(m \times 2)$ matrix defined by
\begin{align*}
H =
\left(\!
\begin{array}{cc}
e^{-k(T_1-t_i)} & 1\\
\phantom{\frac{\frac{1}{1}}{\frac{1}{1}}}\vdots\phantom{\frac{\frac{1}{1}}{\frac{1}{1}}} & \vdots\\
e^{-k(T_m-t_i)} & 1
\end{array}
\!\right),
\end{align*}
and $\pmb{w}_{t_i}$ is an $m$-dimensional normal random vector with zero mean and a constant covariance matrix $V$, to be defined below. Then, based on Equation $(3)$,
\begin{equation}\label{spot_forward_measurement_eqn}
\pmb{y}_{t_i} = \pmb{d}_{t_i} + H \pmb{x}_{t_i} +\pmb{w}_{t_i},
\end{equation}
which is the measurement equation or observation equation. Here $\pmb{w}_{t_i}$ is a $m$-dimensional vector of normal random variables. The $(m \times m)$ covariance matrix $V$ of $\pmb{w}_{t_i}$ is determined, together with the model parameters, as part of the maximum likelihood estimation.
