Using this John Rawnsley: On the universal covering group of the real symplectic group, I found a relatively elegant way to write a partial solution, which makes clear where the sign ambiguity comes from. For this, we use the basis $\hat{\xi}^a=(\hat{q}_1,\dots,\hat{q}_N,\hat{p}_1,\dots,\hat{p}_N)$, such that the ground state covariance matrix $G^{ab}=\langle 0|\hat{\xi}^a\hat{\xi}^b+\hat{\xi}^b\hat{\xi}^a|0\rangle$ is the identity. We then define the associated complex structure and symplectic form
\begin{align}
G\equiv\begin{pmatrix}
1\!\!1 & 0\\
0 & -1\!\!1
\end{pmatrix}\,,\quad J\equiv\begin{pmatrix}
0 & 1\!\!1\\
-1\!\!1 & 0
\end{pmatrix}\,,\quad \Omega\equiv\begin{pmatrix}
0 & 1\!\!1\\
-1\!\!1 & 0
\end{pmatrix}
\end{align}
I use different symbols here, because $J^a{}_b$ is a linear map and $\Omega^{ab}$ is a bilinear form, which only happen to be represented by the same matrix in this specific basis.
The quadratic Hamiltonian is fully characterized by its symmetric coefficient matrix $h_{ab}$, such that $\hat{H}=\frac{1}{2}h_{ab}\hat{\xi}^a\hat{\xi}^b$. With this in hand, we define the symplectic algebra generator $K^a{}_b=\Omega^{ab}h_{bc}$ with the symplectic transformation
\begin{align}
M=e^{t K}=\begin{pmatrix}
A & B\\
C & D
\end{pmatrix}\,.
\end{align}
The expectation value is then given by (I derived it with the method of my previous answer based on the definitions of the Rawnsley paper)
\begin{align}
\langle 0|e^{-i t\hat{H}}|0\rangle=1/\sqrt{\det{\mathcal{C}_M}}\,,
\end{align}
where $\mathcal{C}_M$ is a complex $N$-by-$N$ matrix defined as follows: We first compute $C_M=\frac{1}{2}(M-JMJ)$, which commutes with $J$, i.e., it satisfies $[C_M,J]=0$. However, $C_M$ is a $2N$-by-$2N$ matrix, but because it commutes with $J$ (representing an imaginary unit with $J^2=-1\!\!1$) we can convert it to complex matrix $N$-by-$N$ matrix in the $(+i)$-eigenspace of $J$. We thus find
\begin{align}
C_M=\frac{M-JMJ}{2}=\frac{1}{2}\begin{pmatrix}
A+D & B-C\\
C-B & A+D
\end{pmatrix}\quad\Rightarrow\quad \mathcal{C}_M=\frac{A+D}{2}+i\frac{B-C}{2}\,.
\end{align}
Now we see exactly where the sign ambiguity comes from: The square root function $1/\sqrt{\det{\mathcal{C}_M}}$ has a branch cut and is thus ambiguous up to a sign. Of course, we can always determine the correct sign by hand: We start at $t=0$ with $\langle 0|e^{-i t\hat{H}}|0\rangle=1$ and then increase $t$ to our desired value, while choosing the square root $1/\sqrt{\det{\mathcal{C}_M}}$ in a continuous way.
Example 1. In the special case, where $[K,J]=0$, we also have $[M,J]=0$, such that $C_M=M=e^{Kt}$ with $A=D$ and $B=-C$. Therefore, we have $\mathcal{C}_M=A+i B$. In this case, $K$ is of the following form and we define $\mathcal{K}$
\begin{align}
K=\begin{pmatrix}
0 & k\\
-k & 0
\end{pmatrix}\quad\Rightarrow\quad\mathcal{K}=i k
\end{align}
where $k=\mathrm{diag}(\epsilon_1,\dots,\epsilon_N)$. We then have $\mathcal{C}_M=A+iB=e^{t\mathcal{K}}=e^{itk}$ as complex $N$-by-$N$ matrix. We can then compute the expectation value as
\begin{align}
\langle 0|e^{-i t\hat{H}}|0\rangle=1/\sqrt{\det{\mathcal{C}_M}}=1/\sqrt{\det{e^{it k}}}=e^{-\frac{1}{2}\mathrm{tr}\log{ikt}}=e^{-\frac{it}{2}\mathrm{tr}(k)}=e^{-\frac{it}{2}\sum_i\omega_i}\,.
\end{align}
Here, we got rid of the sign ambiguity by using $\sqrt{\det{e^X}}=e^{\frac{1}{2}\log \det{e^X}}=e^{\frac{1}{2}\mathrm{tr}{\log e^{X}}}$, where we set $\log e^X=X$ without taking branch cuts into account and thereby healing the ambiguity. The result is of course, trivially correct, because for $[J,K]=0$ the state $|0\rangle$ is the ground state of the Hamiltonian with vacuum energy $E_0=\frac{1}{2}\sum_i\omega_i$.
Example 2. The other case that I could treat was the setting where the Hamiltonian $\hat{H}=\sum_i\omega_i(\hat{n}_i+\frac{1}{2})$ has a ground state $|\{0,\dots,0\}\rangle=|0_1\rangle\dots|0_N\rangle$, i.e., it is a tensor product of squeezed 1-mode states, where the initial vacuum $|0\rangle=|0\rangle\dots|0\rangle$ is also a tensor product. Up to a complex phase, we can then relate $|0\rangle$ and $|n_i\rangle$ (arbitrary number excitations of the $n$-th mode of the Hamtiltonian) with a squeezing parameter $r_i$, such that
\begin{align}
|\langle 2n_i|0\rangle|^2=\frac{(2n)! \tanh^{2n_i}{r_i}}{2^{2n_i}(n_i!)^2\cosh{r_i}}\,.
\end{align}
This allows us to use a resolution of the identity, sum the series to find
\begin{align}
\langle 0|e^{-i t\hat{H}}|0\rangle&=\sum_{n_i,n'_j}\langle 0|\{n_1,\dots,n_N\}\rangle\langle\{n_1,\dots,n_N\}|e^{-i \hat{H}}|\{n_1',\dots,n_N'\}\rangle\langle\{n_1',\dots,n_N'\}|0\rangle\\
&=\prod_i\sum_{n_i}|\langle 2n_i|0\rangle|^2e^{-i \omega_i t(2n_i+\frac{1}{2})}\\
&=\prod_i\frac{e^{-\frac{i\omega_it}{2}}\mathrm{sech}(r_i)}{\sqrt{1-e^{-2i\omega_it}\tanh^2{r_i}}}\,.
\end{align}
However, I haven't really tested if the square root in the denominator has the same sign ambiguity.
I have not been able to come up with a smart way to extract the sign in other cases, except manually taking $t$ to its desired value on requiring continuity of $\langle 0|e^{-i t\hat{H}}|0\rangle$ on the way...
Update: I believe that I found a closed form solution of the problem that does not involve integration. It requires the use of the cocycle function $\eta(M_1,M_2)$ defined in John Rawnsley's paper from above (which can be evaluated with elementary matrix operation). In the end, we find that the correct complex phase is given by $e^{i \psi}$ with $\psi=\frac{1}{2}(\mathrm{tr}(uKu^{-1})+\eta(u,e^K)+\eta(ue^{K},u^{-1}))$, where $u$ is a symplectic transformation that brings $K$ into a standard form $\tilde{K}=uKu^{-1}$ in the same basis as $J$ takes its standard form. More precisely, we transform $K$, such that in its real eigenbasis (for complex eigenvalues we keep real blocks with antisymmetric pieces) the complex structure $J$ takes the same form $\Omega$. I may write this up as a more pedagogical note later on, but for now I hope that my solution is useful if somebody has the same question...
