Estimate the Image Using Multi Many Realizations of Its Convolution with a Known Filters Using Wiener Filter Suppose we have a corrupted image $Y = H*X + \epsilon$ that is formed by taking an image $X$, convolving it with a point-spread function $H$, and adding gaussian noise $\epsilon$. Then we know that the Wiener Filter can compute the MMSE estimator of $X$ given $H$ and the signal-to-noise ratio (SNR).
Given a set of $n$ images $Y_i = H_i * X + \epsilon_i$, is there a generalized Wiener Filter estimate of $X$ given the $H_i$'s and the SNRs?

This is a nice question.
The math is actually pretty simple once you embrace the method I derived the Wiener Filter in - How Is the Formula for the Wiener Deconvolution Derived?
So, here is the model:
$$
\boldsymbol{y}_{i} = \boldsymbol{h}_{i} \ast \boldsymbol{x} + \boldsymbol{w}, \; i = 1, 2, \ldots, n
$$
Where $ \boldsymbol{w} $ is an additive white gaussian noise which is independent of the signal and $ \boldsymbol{x} \sim N \left( 0, {\sigma}_{x}^{2} \right) $.
Then the optimization model becomes (The MMSE Estimator):
$$ \arg \min_{\boldsymbol{x}} \frac{1}{ {\sigma}_{n}^{2} } \sum_{i = 1}^{n} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{1}{2 {\sigma}_{x}^{2}} {\left\| \boldsymbol{x} \right\|}_{2}^{2} $$
Where $ H = \begin{bmatrix} {H}_{1} \\ {H}_{2} \\ \vdots \\ {H}_{n} \end{bmatrix} $ is the model matrix where $ {H}_{i} $ is the matrix form of the $ i $ -th filter and $ \boldsymbol{y} = \begin{bmatrix} \boldsymbol{y}_{1} \\ \boldsymbol{y}_{2} \\ \vdots \\ \boldsymbol{y}_{n} \end{bmatrix} $ is a concatenation of the output images in a column stack form.
Then the solution is given by:
$$ \hat{\boldsymbol{x}} = {\left( {H}^{T} H + \frac{ {\sigma}_{n}^{2} }{ {\sigma}_{x}^{2} } I \right)}^{-1} {H}^{T} \boldsymbol{y} $$
