There is a general argument for the existence of the unitary operator. From the Stone-von Neumann theorem, if you know where $a$ is mapped, then there is only one unitary operator that can represent the transformation (up to a global phase).
You can construct this operator more explicitly. One way is to realize that your transformations form a path connected Lie group. Formally, you can view it as the subgroup of matrices of size $2N\times2N$ of the form:
$$
\begin{pmatrix}
A & B\\
B^* & A^*
\end{pmatrix}
$$
with your additional constraints:
$$
A A^\dagger -BB^\dagger =1\\
AB^T-BA^T =0
$$
You can consider a path from the identity to your transformation, namely $A_t,B_t$ with $t\in[0,1]$ and:
$$
\begin{align}
A_0&=1 & B_0&=0 \\
A_1&=A& B_1&=B
\end{align}
$$
This gives a varying $a_t$ satisfying the ODE:
$$
\dot a_t= \alpha_ta_t+ \beta_ta_t^\dagger
$$
with the matrices $\alpha_t,\beta_t$ defined by:
$$
\frac{d}{dt} \begin{pmatrix}
A_t & B _t\\
B _t ^* & A _t ^*
\end{pmatrix} = \begin{pmatrix}
\alpha_t& \beta_t\\
\beta_t^* & \alpha_t^*
\end{pmatrix} \begin{pmatrix}
A_t & B _t\\
B _t ^* & A _t ^*
\end{pmatrix}
$$
You want to recognize the ODE as a Heisenberg equation from a suitably chosen quadratic Hamiltonian:
$$
\dot a_t =-i[a,H]
$$
Using the CCR’s you can check that a suitable candidate (unique up to additional multiple of identity) is:
$$
H_t= i(\alpha_t)_{ij}a_i^\dagger a_j+\frac{i}{2}(\beta_t)_{ij}a_i^\dagger a_j^\dagger-\frac{i}{2}(\beta_t^*)_{ij}a_i a_j
$$
$H_t$ is hermitian and gives the correct equations of motion using the relations (taking the derivatives of the constraints):
$$
\alpha_t+\alpha_t^\dagger =0 \\
\beta_t-\beta_t^T=0
$$
Your unitary operator is given by the time ordered exponential:
$$
U=\mathcal T \exp\left(\int_0^1dt H_t\right)
$$
Actually, it turns out that the exponential is surjective on the group, so you can choose
Choose a path with constant $\alpha_t,\beta_t$. This gives a constant $H_t$ and the time ordered exponential is a simple exponential.
Physically, you are constructing the unitary operator from harmonic oscillators and squeeze operators. The process is explicit in the case $N=1$.
Hope this helps.
