The following does not fully resolve your question, but does supply an unbiased estimator for the mean, based on an approximation to the optimal aggregation of samples. This may prove useful for analysing certain properties of the estimator, both in theory and practice, particularly in the limit of large sample size.
For the $K$ distributions, write their laws as $\{ \pi_k \}_{k = 1}^K$, and write
\begin{align}
\mathbf{E}_{\pi_k} \left[ X \right] &= \mu \\
\textbf{Var}_{\pi_k} \left[ X \right] &= \sigma^2_k
\end{align}
Suppose now that for each $k$, you draw $N_k \geqslant 3$ samples from $\pi_k$, and label them as $\{ X_k^i \}_{i = 1}^{N_k}$.
A short calculation analogous to your own indicates that if $\{ \sigma_k^2 \}_{k = 1}^K$ were all known, then the optimal way of aggregating all of the samples into an unbiased estimator would be
\begin{align}
\alpha_k &\triangleq \frac{N_k / \sigma_k^2}{\sum_{\ell = 1}^K N_\ell / \sigma_\ell^2}\\
\hat{\mu}_k &\triangleq \frac{1}{N_k}\sum_{k = 1}^K X_k^i \\
\hat{\mu} &\triangleq \sum_{k = 1}^K \alpha_k \hat{\mu}_k \\
&= \sum_{k = 1}^K \sum_{i = 1}^{N_k} \frac{\alpha_k}{N_k} X_k^i
\end{align}
Now, the natural thing is to attempt to estimate the $\alpha_k$ using your collection of samples. The "obvious" solution here (replace the $\sigma_k^2$ by estimates built from the samples) will in general lead to a biased estimator, as the estimate of $\alpha_k$ will depend on $X_k^i$, and so the expectation of their product will be more complicated to control.
A resolution to this is: for each pair $(k, i)$ such that $1 \leqslant i \leqslant N_k$, form an estimate $\hat{\alpha}_k^i$ of $\alpha_k$ which is independent of $X_k^i$ .
One such approach would be the following:
- For $1 \leqslant k \leqslant K$, form the following estimate of $\sigma_k^2$:
$$\hat{\sigma}_k^2 \triangleq \frac{1}{N_k - 1} \sum_{i = 1}^{N_k} \left( X_k^i - \hat{\mu}_k \right)^2.$$
This estimator happens to be unbiased for $\sigma_k^2$, but this fact is not used later on. Other choices may be more appropriate.
- For $1 \leqslant k \leqslant K, 1 \leqslant i \leqslant N_k$, form the following "leave-one-out" estimate of $\sigma_k^2$ as follows:
\begin{align}
\hat{\mu}_{k, -i} &\triangleq \frac{1}{N_k - 1} \sum_{j \neq i} X_k^j\\
\hat{\sigma}_{k, -i}^2 &\triangleq \frac{1}{N_k - 2} \sum_{j \neq i} \left( X_k^i - \hat{\mu}_{k, -i} \right)^2.
\end{align}
Note that these can be computed quickly by noting that
\begin{align}
\hat{\mu}_{k, -i} &= \hat{\mu}_k + \frac{\hat{\mu}_k - X_k^i}{N_k - 1} \\
\hat{\sigma}_{k, -i}^2 &= \hat{\sigma}^2_k + \frac{\hat{\sigma}^2_k - \left(\hat{\mu}_k - X_k^i\right)^2}{N_k - 2}
\end{align}
Again, this estimator happens to be unbiased for $\sigma_k^2$, but this fact is not used later on, and other choices may be more appropriate. Crucially, $\hat{\sigma}_{k, -i}^2$ is independent of $X_k^i$.
- For $1 \leqslant k \leqslant K, 1 \leqslant i \leqslant N_k$, form the following "local" estimate of $\alpha_k$
$$\hat{\alpha}_k^i \triangleq \frac{N_k / \hat{\sigma}_{k, -i}^2 }{\sum_{\ell \neq k} N_\ell / \hat{\sigma}_\ell^2 + N_k / \hat{\sigma}_{k, -i}^2},$$
which is independent of $X_k^i$.
One could then report the final unbiased estimate of $\mu$ as
$$\hat{\mu} = \sum_{k = 1}^K \sum_{i = 1}^{N_k} \frac{\hat{\alpha}_k^i}{N_k} X_k^i,$$
which will be unbiased, as advertised.
At least asymptotically, one should be able to say that the $\hat{\alpha}_k^i$ will be close to the true $\alpha_k$. The reasoning here is roughly that $\hat{\sigma}_{k, -i}^2, \hat{\sigma}^2_k$ each converge to $\sigma_k^2$ (in probability, almost surely, etc.), and so the continuous mapping theorem should apply.
With this in mind, the variance should asymptotically be the same as if we knew the values of $\{ \sigma_k \}_{k = 1}^K$ a priori, and it is not too difficult to show that this must dominate the naive strategy.