## 3 Materials and Methods
### Recursive Estimation and Parameters Settings  
Our algorithm has a lot of prior variables, as shown in the previous derivation; mainly \(\mathbf{K},\mathbf{k}(x),\mathbf{P}_{0},\mathbf{R}_{t}\), and \(\mathbf{F}\). First, the parameter selection of \(\mathbf{R}_{t}\) is mainly based on the measurement error of the observation itself. Therefore, the observed quantities are independent of each other. Therefore, for \(\mathbf{R}_{t}\), we need the corresponding probability distribution to obtain the variance.  
Second, we can find the value of \(\mathbf{P}_{0}=T\mathbf{K}\) using the length of the sample data, or we can just specify the value of \(\mathbf{P}_{0}=0\) or a smaller value, when \(T\) is large.  
Third, we must specify drift functions to obtain \(\mathbf{f}(x)\) and \(\mathbf{F}\), mainly to reflect spatial non-stationarity. However, as the specific order of the spatial trend distribution cannot be obtained, the lower order is generally chosen to reflect the possible spatial trend characteristics, and the experiment proves that this choice also obtains good results. In fact, this phenomenon also conforms to the characteristics of universal Kriging estimation. In the process of universal Kriging estimation, it is necessary to calculate the variogram or covariance function of the residual error by removing the drift term. Compared with the calculation without drift term, the result can be made more accurate by removing the trend term of lower order.  
Finally, it is difficult to choose the general \(k(x_{a},x_{b})\) in Universal Kriging theory, to determine \(\mathbf{K}\) and \(\mathbf{k}(x)\). To simplify the calculation, we chose intrinsic random functions (IRF), which are a more generalized covariance function. In this paper, we propose two generalized covariance functions. [PERSON] [15] proposed an IRF model of order \(k\), with expression \(k(\parallel x_{a}-x_{b}\parallel_{2})=\parallel x_{a}-x_{b}\parallel_{2}^{2k+1}\). [PERSON] [16] proposed an IRF function similar to a spline, with expression \(k(\parallel x_{a}-x_{b}\parallel_{2})=\parallel x_{a}-x_{b}\parallel_{2}^{2} \log\parallel x_{a}-x_{b}\parallel_{2}\). Therefore, the form of \(k(x_{a},x_{b})\) in this paper mainly refers to these two forms, and good results were obtained in the experimental verification. The specific steps are provided in Algorithm 1. As shown in table, this FUKSS model maintains a fixed-size state variable \(c_{t}\) and transfer matrix \(P_{t}\). The complexity of calculation at time \(t\) is only related to the quantity of valid sampling data \(Y_{t}\).  
```
Initialize: \(t=1,P_{0}=\alpha K\), \(\epsilon_{0}=\mathbf{0}(n*1)\).  Given: \(\mathbf{K}\), \(\mathbf{k}(x)\), \(\mathbf{\alpha}\), \(\mathbf{R}_{t}\), \(\mathbf{F}\), \(\mathbf{W}_{t}\)  State equations: \[Z_{t}(x)=Z_{t-1}(x)+f^{T}\mathbf{\beta}_{t}+\varepsilon_{t}(x)\]  Observation equation: \[Y_{t}(x)=Z_{t}(x)\] For \(t=1,\dots,T\)  Calculate the matrices: \[\mathbf{L}_{t}=(\mathbf{W}_{t}\mathbf{P}_{t-1}\mathbf{W}_{t}^{T}+\mathbf{W}_{t}\mathbf{K} \mathbf{W}_{t}^{T}+\mathbf{R}_{t})^{-1}\] \[\mathbf{M}_{t}=(\mathbf{F}\mathbf{W}_{t}^{T}\mathbf{L}_{t}\mathbf{W}_{t}\mathbf{F}^{T})^{ -1}\mathbf{F}\mathbf{W}_{t}^{T}\mathbf{L}_{t}\] \[\mathbf{G}_{t}=(\mathbf{P}_{t-1}+\mathbf{K})\mathbf{W}_{t}^{T}\mathbf{L}_{t}(\mathbf{I}- \mathbf{W}_{t}\mathbf{F}^{T}\mathbf{M}_{t})\]  Calculate the desired estimate using: \[\widehat{Y}_{t}(x)=\mathbf{k}^{T}(x)\mathbf{K}^{-1}\widehat{\varepsilon}_{t-1}+(\mathbf{k }^{T}(x)\mathbf{K}^{-1}\mathbf{G}_{t}+\mathbf{f}^{T}(x)\mathbf{M}_{t})(\mathbf{Y}_{t}-\mathbf{W}_{t} \widehat{\varepsilon}_{t-1})\] Updated parameters : \[\mathbf{P}_{t}=(\mathbf{P}_{t-1}+\mathbf{K})(\mathbf{I}-\mathbf{W}_{t}\mathbf{G}_{t}^{T})\] \[\widehat{\varepsilon}_{t}=\widehat{\varepsilon}_{t-1}+\mathbf{G}_{t}( \mathbf{Y}_{t}-\mathbf{W}_{t}\widehat{\varepsilon}_{t-1})\] end
```  
**Algorithm 1**FUKSS Algorithm