I set $\hbar = 1$
Another way to come at this is to consider the Hamiltonian to be the generator of unitary time evolution. That is, say we have an initial state $|\psi(0)\rangle$. We can take it as a postulate of quantum mechanics that the state at a later time is given by
$$
|\psi(t)\rangle = U(t) |\psi(0)\rangle
$$
Where $U$ is a unitary operator. By taking the time derivative of this expression we arrive at the Schrodinger equation
$$
\frac{d}{dt}|\psi(t)\rangle = -i H(t) |\psi(t)\rangle
$$
Where I've defined the Hamiltonian
$$
H(t) = i \frac{dU}{dt} U^{\dagger}
$$
From this definition and unitary of $U(t)$ it can be proven that
$$
H(t) = H^{\dagger}(t)
$$
and, if $H(t)$ happens to be time-independent:
$$
U(t) = e^{-i H t}
$$
What you seem to be gesturing at in your question is that the Hermiticity of $H$ seems to intuitively follow from the reversibility of dynamics/interactions. This reversibility is exactly encoded in the unitarity of the time evolution operator, so deriving Hermiticity of $H$ from this unitarity should then, maybe, be satisfying to you.
I might as well collect my notational simplification in a non-answer, instead of in garlands of comments...
Define $\sqrt{G}\equiv \cosh g \leadsto \sqrt{G-1}= \sinh g$. Moreover, focus on the states $|n\rangle\equiv |n,n\rangle$. Your "hamiltonian" keeps you in this class, and notation convention,
$$
(b^\dagger b - \sinh^2 g)|n\rangle \\ = (\cosh 2g)~ n |n\rangle + \tfrac{1}{2} \sinh 2g \bigl ((n+1)|n+1\rangle +n|n-1\rangle\bigr ) ,
$$
which you may write as the standard Fock infinite-dimensional matrix like this and contemplate the eigenvectors of,
$$
M= \begin{pmatrix}
0& \tfrac{1}{2} \sinh 2g & 0 & ...\\ \tfrac{1}{2} \sinh 2g&
\cosh 2g & \sinh 2g &0& ... \\
0& \sinh 2g & 2 \cosh 2g & \tfrac{3}{2} \sinh 2g & 0& ... \\ 0&0&...
\end{pmatrix}.
$$
Edit in response to question's 3rd edit:
To find the eigenvectors of the matrix M, you write the eigenvalue conditions, which specify a recursive determination of their coefficients, given the escalated structure of the Matrix, and so the (infinity) of eigenvectors will all be infinite-dimensional. (But, unlike the simple eigenvectors of $\hat X$ in Fock space, $ \frac{1}{\pi^{1/4}} e^{-x^2/2} e^{\sqrt{2} xa^\dagger} e^{-a^{\dagger ~2}/2} =\frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle
$, they may be very complicated...)
The recursive determination of the eigenvectors goes as follows,
$$
M|\psi\rangle= \lambda |\psi\rangle, ~~~ |\psi\rangle = \sum_n c_n|n\rangle ~~\leadsto \\
\lambda c_0= \tfrac{1}{2} \sinh 2g ~~c_1, \\
\lambda c_1 = \tfrac{1}{2} \sinh 2g ~~c_0 +
\cosh 2g ~~ c_1 + \sinh 2g ~~ c_2, \\ ... ,
$$
so you may recursively solve for $c_1$ in terms of $\lambda, ~g$ and $c_0$, then for $c_2$ in terms of the above, etc... Conceptually, the infinity of $c_n$ are uniquely determinable in terms of $\lambda, ~g$ and $c_0$, given the stepwise nature of the matrix M ! If the answer were a simple closed form, you could advance to issues of convergence, etc..., as in the linked example., but only then: don't let such qualms stop you.
For a zero eigenvalue, taking $c_0=1$ to work the normalization later, you find $(1,0,-1/2, 2/3 ~ \coth 2g , 3/8-\coth^2 2g,...)^T$, which is to say
$$
|\psi\rangle=
|0\rangle -\tfrac{1}{2}|2\rangle + \tfrac{2}{3} ~ \coth 2g |3\rangle +(\tfrac{3}{8}-\coth^2 2g)|4\rangle +...
$$
Best Answer
$\vert\omega_1\rangle$ would not be a basis, but the set $\{\vert\omega_1\rangle,\vert\omega_2\rangle, \ldots,\vert\omega_n\rangle\}$ is a basis. The eigenvector $\vert\omega_i\rangle$ is such that $\Omega\vert\omega_i\rangle=\omega_i\vert\omega_i\rangle$ so...
With the identification $$ \vert\omega_1\rangle \to \left(\begin{array}{c} 1 \\ 0 \\ \vdots \\ 0\end{array}\right)\, ,\quad \vert\omega_2\rangle \to \left(\begin{array}{c} 0 \\ 1 \\ \vdots \\ 0\end{array}\right)\, ,\ldots \, , \vert\omega_n\rangle \to \left(\begin{array}{c} 0 \\ \vdots \\ 0 \\ 1 \end{array}\right) $$ and the matrix representation of $\Omega$ as you suggest you find by simple matrix multiplication that $$ \Omega \vert \omega_n\rangle = \omega_n\vert\omega_n\rangle $$ as per the properties of eigenstates of $\Omega$. Note that under this identification the vectors $\vert\omega_i\rangle$ are an orthonormal basis since $\langle \omega_i\vert\omega_j\rangle=\delta_{ij}$.