Here is a calculation based on a general Bogoluibov transformation
$$
a_i= u_{i\beta}b_\beta +v^*_{i\beta}b^\dagger_\beta\\
a^\dagger_i= v_{i\alpha} b_\alpha +u^*_{i\alpha}b^\dagger_\alpha.
$$
With $\vert 0 \rangle_b$ being the vacuum state for the $b_\alpha$'s we have
$$
N=\sum_i { _b\langle{0}\vert{ a^\dagger_i a_i}\vert{0}\rangle _b} \\
={\sum_{i,\alpha,\beta} { _b\langle{0}\vert{( v_{i\alpha} b_\alpha +u^*_{i\alpha}b^\dagger_\alpha)(u_{i\beta}b_\beta +v^*_{i\beta}b^\dagger_\beta)}\vert{0}\rangle_b}}\\
= \sum_{i,\alpha,\beta} (v_{i\alpha} v^*_{i\beta}) {_b\langle{0}\vert{ b_\alpha b^\dagger_\beta }\vert{0}\rangle_b}\\
=\sum_{\alpha,\beta} v_{i\alpha} v^*_{i\beta} \delta_{\alpha\beta}\\
= \sum_\alpha |v_\alpha|^2.
$$
We can also write
$$
N= \frac 12 \sum_{\alpha=1}^N\left((|v_\alpha|^2-|u_\beta|^2) + (|u_\alpha|^2+|v_\alpha|^2)\right)\\
=\frac 12 \sum_{\alpha=1}^N\left((|v_\alpha|^2-|u_\beta|^2) + 1\right),
$$
which looks like the contribution of the Dirac sea negative-energy Bogoliubov-de Gennes equation eigenstates $(v^*_\alpha, u^*_\alpha)$, but corrected by an infinite normal-ordering counterterm (the sum of all the ``+1''s).
As an aside, note I'm not sure that your BCS ground state is normalized. I like to write
$$
\vert{0}\rangle _b ={\mathcal N} \exp\left\{\frac 12 a^\dagger_i a^\dagger_jS_{ij}\right\} \vert{0}\rangle_a
$$
where
$$
S_{ij}= v^*_{i\alpha}(u^*)^{-1}_{\alpha j}
$$
is skew symmetric
and
$$
{\mathcal N}= {\rm det}^{-1/4}(I+S^\dagger S)
$$
is derived by using the fermionic version of MacMahon's master theorem. This way makes clear that the BCS ground state is a coherent state superposition of Cooper pairs with pair "wavefunction" $S_{ij}$.
Note added: Yes your state is normalized. It's just not the way I like to write it.
