How to initialize observation Matrix in Kalman Filter when there is no clear relationship between measurement and state?
In order to get a better understanding I used the method proposed in the paper to build the measurement function $y=h(x)$ using a Matrix $H$ and the state transition function $x_{n+1}=f(x_n)$ using a Matrix $A$. I did this for a practical example to get a feeling how it works.
To be clear in the paper you cited they call the measurements $z$ instead of $y$. So their measurement function is:

Their state transfer function is:

If the physical relationship between the measurement $z$ and the state $x$ is unknown one can use the formulas from the paper to calculate $H$ and $A$ depending on past states and measrements.
The example I tried is from Plane with constant speed. A radar system measures the distance to a flying plane that travels with constant speed. Starting position is 30 km away. Starting speed is 40 m/s. The state Vector has two values: position and speed of the airplane. So the inital state would be: $x_1=[30000, 40]$. Because the radar only measures the distance, the measurement equation simply is: $z_n = \begin{pmatrix} x_{1,n} \\ 0 \end{pmatrix}$. Thats because only the distance ($x_{n,1}$) gets measured but not the velocity ($x_{n,2}$).
The state transition equation is $\begin{pmatrix} x_{1,n+1} \\ x_{2,n+1}\end{pmatrix} = \begin{pmatrix} x_{1,n} + \Delta t \cdot x_{2,n} \\x_{2,n}\end{pmatrix}$. The velocity $x_{n,2}$ stays the same between step $n$ and $n+1$ but the distance changes by the multiplication of the time step $\Delta t$ and the speed during that time step $x_{2,n}$.
Using those equations and the measurement values provided in the example I get the following:

Next we try to achieve the same without writing the functions $f$ and $h$ explicitly but we try to use the formulas from the paper you cited.
For that also the recorded states $x_{n}$ and $x_{n+1}$ are needed. Luckily they are also provided in the example:

Using the Formulas for $H$ and $A$ and the recorded traces it is possible to calculate them.

$H$ in the case of our example becomes a $1\times2$ Matrix. Thats because $z$ is a $1\times10$ vector and $x_n$ is a $2\times10$ vector. As the latter is transposed $H$ ends up being $1\times2$.
Quite the same way $A$ ends up being a $2\times2$ matrix. $x_{n+1}$ and $x_n$ are both a $2\times10$ matrix. As the latter is transposed by the multiplication one ends up with a $2\times 2$ matrix.
Now we can rewrite the equations $f$ and $h$.
$$x_{n+1}=f(x_n): x_{n+1}=A\cdot x_n$$
and
$$z_n=h(x_n): z_n=H\cdot x_n$$
If I put in the recorded values I get $H=\begin{pmatrix} 0.956 &amp; 35.38 \end{pmatrix}$ and $A=\begin{pmatrix} 1.00 &amp; 4.98 \\ 0.00 &amp; 1.00 \end{pmatrix}$.
With that I rerun the Kalman filter and I get the following:

This looks very similar to the former result which was derived using a physical modell. So in this case the method from the paper seems to work quite fine.
I did this using Matlab. So in case you are interested in the code:
close all

initialStateGuess = [30000 40];
n=length(initialStateGuess);
dt=5;

x_state=[30182  30351.4 30573.3 30769.5 31001.5 31176.4 31333.2 31529.4 31764.3 31952.9;38.2    36  40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
x_state2=[30373 30531.6 30774.3 30968.1 31216.8 31371.5 31509.2 31715.4 31974.8 32152.4; 38.2   36.0    40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
yMeas = [30110 30265 30740 30750 31135 31015 31180 31610 31960 31865];% Measurement Values
H=(yMeas*x_state2')/(x_state2*x_state2');
A=(x_state2*x_state')/(x_state*x_state');

% h=@(x)[x(1)];% Measurement function
% f=@(x)[x(1)+dt*x(2);x(2)];% State transition function

h=@(x)[H(1)*x(1)+H(2)*x(2)];
f=@(x)[A(1,1)*x(1)+A(1,2)*x(2);A(2,1)*x(1)+A(2,2)*x(2)];

ukf = unscentedKalmanFilter(...
    f,... % State transition function
    h,... % Measurement function
    initialStateGuess,...
    'HasAdditiveMeasurementNoise',true,...
    'HasAdditiveProcessNoise',true);

yTrue=linspace(30000,32000,10);
R = var(yMeas-yTrue);%0.2; % Variance of the measurement noise v[k]

Q=var(x_state-yTrue);
ukf.MeasurementNoise = R;
ukf.ProcessNoise = 30;

Nsteps = length(yMeas); % Number of time steps
xCorrectedUKF = zeros(Nsteps,n); % Corrected state estimates
PCorrected = zeros(Nsteps,n,n); % Corrected state estimation error covariances
e = zeros(Nsteps,1); % Residuals (or innovations)

for k=1:Nsteps
    % Let k denote the current time.
    %
    % Residuals (or innovations): Measured output - Predicted output
    e(k) = yMeas(k) - h(ukf.State); % ukf.State is x[k|k-1] at this point
    [xCorrectedUKF(k,:), PCorrected(k,:,:)] = correct(ukf,yMeas(k));
    predict(ukf);
end

plot([5:5:50],xCorrectedUKF(:,1))
hold on
plot([5:5:50],yMeas)
plot([5:5:50],yTrue)
legend('Kalman FIlter state','Measurement Values','True Value')
grid on