Mixed models and geostatistics : interpretation

Hi, I thank in advance all persons who can help me to understand my problem more clearly.

let's say we have 20 air quality sensors that crisscross a city. Resulting in N observations (each sensor has n_k observations), each observation z is linked to the coordinates where it was collected.

first, I use a gaussian process to model my data: I assume that the hidden process Y (which represents the dispersion of the pollutants) follows a gaussian process with an exponential function as covariance function, and that the observations of the sensors are just noisy independent observations, classical in geostatistics :

Z/Y ~ N( Y, σI )

Y ~ N( µ , R( ))

Where µ can be a constant or a regression with other covariate collected by the sensors ( ex : proximity to highway ) ,and R() is an exponential function of the distance depending on several parameters , and σ represent the mesurment error. Resulting in the marginal model :

Z ~ N( µ , R( ) + σI)

This model work just fine, but does not take into account the fact that my data comes from unreliable ,different sensors.

I want my model to be able to have several intercepts in Z for each sensor k, this is where i came across linear mixed models :

Z_k/Y ~ N( Y + a_k, σI )

Y ~ N( µ , R( ))

a_k ~N( 0 , γ )

I understand that the a_k must be centred at 0, or they will not be distinguished from µ , and the model will be unidentifiable.

In the first model (geostatistics), the resulting covariance matrix Σ is equal to R(xi,xj) for any input not in the diagonal, and R(xi,xi) + σ in the diagonal.

Integrating out of the a_k in the second model results in a more complex covariance matrix Σ , and we do not directly observe the a_k.

Here Σ is equal to R(xi,xj) for every 2 observations not from the same sensor, R(xi,xj) + γ for every two observations from the same sensor and R(xi,xi) + σ + γ on the diagonal. (is this correct ?)

My understading is that γ capture all the variation within the sensors.

My questions are : 1) Is my second model correct?

2 ) We have seen that the γ present between 2 observations from the same sensor in the covariance matrix, can be interpreted as different intercepts in the model for each sensor k, centered on 0. What is the interpretation if, instead of putting the same parameter, I put γ_k in the covariance matrix for each sensor k, thus adding (k-1) parameter (each sensor has its own variance).

3) Same question if i simply put different σ_k on the diagonal for each observation from sensor k.

If anyone has a lead or an idea, a paper or a book chapter that could help me, I would be very grateful.

Thank you.