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}{ 2 {\sigma}_{n}^{2} } \sum_{i = 1}^{n} {\left\| {H}_{i} \boldsymbol{x} - \boldsymbol{y}_{i} \right\|}_{2}^{2} + \frac{1}{2 {\sigma}_{x}^{2}} {\left\| \boldsymbol{x} \right\|}_{2}^{2} = \arg \min_{\boldsymbol{x}} \frac{1}{ { 2 \sigma}_{n}^{2} } {\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} $$
