Multivariate linear regression with dependant observations the Multivariate linear regression model is given by
$\underset{n \times d}{\mathbf{Y}} = \underset{n \times k}{\mathbf{X}} \hspace{2mm}\underset{k \times d}{\mathbf{B}} + \underset{n \times d}{\mathbf{E}}$ 
With $\mathbf{Y}$ is data, $\mathbf{X}$ is a non stochastic design matrix, $\mathbf{E}$ are normally distributed model errors, and $\mathbf{B}$ are the regression coefficients to be estimated. Also $n$ is the number of observation (e.g.,time samples in my specific context), $d$ the number of variables (e.g., time series at different pixels of an image) and $k$ the number of explanatory variables.  
All the books i read until now treat this problems using the assumption that the obervations (the different rows of $\mathbf{E}$) are independants whereas the different variables (the different columns of $\mathbf{E}$) are dependant with the $d \times d$ covariance structure $\boldsymbol\Sigma_d$.
In this context, the model can be reparametrized as a multiple linear regression model:
$\underset{nd \times 1}{\mathbf{y}} = \underset{nd \times nk}{\mathbf{D}}  \hspace{2mm} \underset{nk \times 1}{\boldsymbol{\beta}} + \underset{nd \times 1}{\mathbf{e}}$ with $\mathbf{e}\sim \mathcal{N}(\mathbf{0},\mathbf{I}_n \otimes\boldsymbol\Sigma_d$) where $\otimes$ is the kronecker product and solved using Generalized Least Square or Feasible GLS if $\boldsymbol\Sigma_d$ is unknown. 
QUESTIONS: 
1-What happens when the observation are not independants ?. 
2-Can we reparametrize the model using  with $\mathbf{e}\sim \mathcal{N}(\mathbf{0},\mathbf{\Omega}_n \otimes\boldsymbol\Sigma_d$)  and solve it using the GLS when $\mathbf{\Omega}_n$ and $\boldsymbol\Sigma_d$ are known .
3- Why this general problem  (non independance of observations) is not treated in textbooks . Is it because the FGLS is too complex (i.e., find an estimate of both  $\mathbf{\Omega}_n$ and $\boldsymbol\Sigma_d$)
Best
 A: 1. If the rows of $\mathbf{Y}$ are not independent, there is in general nothing one can say about the structure of $\mathbf{y}:=\mathrm{vec}(\mathbf{Y})$ other than $\mathbf{y}\sim N\left((\mathbf{I} \otimes \mathbf{X})\mathrm{vec}(\mathbf{B}), \mathbf{\Psi}\right)$ for some covariance matrix $\mathbf{\Psi}$, which is the classical GLS model.
2. Yes you can, given that you are willing to assume this more restrictive structure and that you also know the matrices. Regarding restrictions, you may notice that this decomposition implies that the autocovariance function for any one pixel is proportional to that of any other.
If you are indeed willing to make this assumption and know the matrices, then just set $\mathbf{\Psi} = \mathbf{\Omega}_n \otimes \mathbf{\Sigma}_d$ and your GLS estimator of $\mathbf{b} = \mathrm{vec}(\mathbf{B})$ is given by $\mathbf{Z}(\mathbf{Z}'\mathbf{\Psi}^{-1}\mathbf{Z})^{-1}\mathbf{Z}'\mathbf{\Psi}^{-1}\mathbf{y}$, where $\mathbf{Z} = \mathbf{I} \otimes \mathbf{X}$.
3. I believe your guess at this point is pretty much correct. The model you suggest for $\mathbf{e}$ is equivalent to assuming a matrix normal distribution for $\mathbf{E}$. This paper considers estimation of the component covariance matrices $\mathbf{\Omega}_n$ and $\mathbf{\Sigma}_d$ using maximum likelihood. Paywall link to published version in Journal of Multivariate Analysis, 2016. 
Update:
There is in fact a newer paper (here is an arxiv version) which settles the uniqueness and existence question for maximum likelihood in this model.
Their Theorem 1 in section 4 gives the following result. Suppose we observe $N$ independent copies of the matrix normal variable $\mathbf{Y}$ and $\mathbf{X} = \mathbf{1}_n$. Then:


*

*If $N < \max(n/d, d/n) + 1$ the MLE of $\mathbf{\Psi}$ does not exist.

*If $N > d/n + n/d + 1$ the MLE exists and is unique almost surely.

*For $N \in [\max(n/d, d/n) + 1, d/n + n/d + 1]$ there are examples of when the MLE exists and when it does not, so no general statement can be made.


In particular, note that your case corresponds to having just $N=1$ observation.
