## 2 Related Works
### The Updated Algorithm Considering Vertical Turbulent Wind Field  
The essence of the non-negative matrix factorization (NMF) is a bound-constrained minimization problem. As for a non-negative matrix \(\mathbf{L}\in R^{M\times N}\) and a pre-specified positive integer \(P<\min\{M,N\}\), NMF can be utilized to find two non-negative matrices \(\mathbf{W}\in R^{M\times P}\) and \(\mathbf{V}\in R^{P\times N}\) satisfying [28]:  
\[\mathbf{L}\approx\mathbf{W}\mathbf{V} \tag{18}\]  
where \(\mathbf{W}\mathbf{V}\) is the NMF of \(\mathbf{L}\). The conventional approach proposed by [PERSON] et al. [29] can be utilized to find the NMF of \(\mathbf{L}\):  
\[\min_{W,V}f(W,V)\equiv\frac{1}{2}\sum_{i=1}^{M}\sum_{j=1}^{N}\left(L_{ij}-(WV) _{ij}\right)^{2},\;W_{im}\geq 0,\;V_{nj}\geq 0,\;\forall i,m,n,j \tag{19}\]  
It is noted that Equation (19) is a standard optimization problem with bound constraint. The simple iteration algorithm proposed by [PERSON] and [PERSON] [30] is the widely-used approach to solve Equation (19). \[\mathbf{W}_{ip}\leftarrow\mathbf{W}_{ip}\frac{\left(\mathbf{S}\mathbf{V}^{T} \right)_{ip}}{\left(\mathbf{W}\mathbf{V}\mathbf{V}^{T}\right)_{ip}}\quad \mathbf{V}_{pj}\leftarrow\mathbf{V}_{pj}\frac{\left(\mathbf{W}^{T}\mathbf{S} \right)_{pj}}{\left(\mathbf{W}^{T}\mathbf{W}\mathbf{V}\right)_{pj}} \tag{20}\]  
At every iteration, \(W_{ip}\) and \(V_{pj}\) are multiplied by certain factors. When the iteration is stopped, \(\mathbf{W}\) and \(\mathbf{V}\) are strictly positive. It is noted that Equation (19) will become invariant if the stationary points of the \(\mathbf{W}\) and \(\mathbf{V}\) are reached in Euclidean distance [30]. The \(H^{jm}(\omega,t)\) can be further factorized using the NMF, then the FFT technique can be utilized in the non-stationary SRM. \[H^{jm}_{M\times N}(\omega,t)\approx\sum_{q=1}^{P}w_{q}^{jm}(t)v_{q}^{jm}(\omega) \tag{21}\]Substituting Equation (21) into Equation (13) and using [PERSON]'s formula satisfying:  
\[x_{j}(t)=\text{Re}\Bigg{\{}2\sum_{m=1}^{j}\sum_{q=1}^{P}w_{q}^{jm}(t)\sum_{l=1}^{N }v_{q}^{jm}\big{(}\omega_{l}\big{)}\sqrt{\Delta\omega}e^{i(\omega_{l}t-\delta_{jm }(\omega_{l}t)+\Phi_{ml})}\Bigg{\}} \tag{22}\]  
Note that the simulation efficiency of Equation (22) is much higher when the coherence matrix \(\Gamma(\omega)\) is time-invariant. From Equation (3), \(\mathbf{S}(\omega,t)\) can be rewritten as [31]:  
\[\mathbf{S}(\omega,t)=\mathbf{D}(\omega,t)\mathbf{\Gamma}(\omega)\mathbf{D}^{T} (\omega,t) \tag{23}\]  
\[\mathbf{D}(\omega,t)=\text{diag}\Bigg{[}\sqrt{\mathcal{S}_{11}(\omega,t)}, \sqrt{\mathcal{S}_{22}(\omega,t)},\cdots,\sqrt{\mathcal{S}_{nn}(\omega,t)} \Bigg{]} \tag{24}\]  
\[\mathbf{\Gamma}(\omega)=\left[ \begin{array}{cccc}1&\gamma_{12}(\omega)& \cdots&\gamma_{1n}(\omega)\\ \gamma_{21}(\omega)&1&\cdots&\gamma_{2n}(\omega)\\ \vdots&\vdots&\ddots&\vdots\\ \gamma_{n1}(\omega)&\gamma_{n2}(\omega)&\cdots&1\end{array} \right] \tag{25}\]  
It is obvious that \(\Gamma(\omega)\) is a Hermitian matrix, which also has the positive property. When the elements in \(\Gamma(\omega)\) are complex, the EPSD matrix can be rewritten in polar form [32]:  
\[\mathbf{S}(\omega,t)=\big{|}\mathbf{S}(\omega,t)\big{|}\otimes\mathbf{\Theta }(\omega) \tag{26}\]  
\[\big{|}\mathbf{S}(\omega,t)\big{|}=\left[ \begin{array}{cccc}S_{11}(\omega,t )&\big{|}\mathcal{S}_{12}(\omega,t)\big{|}&\cdots&\Big{|}\mathcal{S}_{1n}( \omega,t)\\ \mathcal{S}_{21}(\omega,t)&S_{22}(\omega,t)&\cdots&\Big{|}\mathcal{S}_{2n}( \omega,t)\Big{|}\\ \vdots&\vdots&\ddots&\vdots\\ \mathcal{S}_{n1}(\omega,t)&\big{|}\mathcal{S}_{n2}(\omega,t)\big{|}&\cdots&S_ {nn}(\omega,t)\end{array} \right] \tag{27}\]  
Note that the phase angle in Equation (12) was the function coupled with time and frequency. However, it is only determined by frequency with regard to Equation (23). Thus, \(\Theta\) can be expressed as:  
\[\mathbf{\Theta}(\omega)=\left[