I am not sure that you want more formal explanation but let me sketch. To develop GL expansion by order parameter $\Delta$, one should start from original 4-fermion theory with contact attraction (BCS theory). Then one should perform decoupling in Cooper channel which means Hubbard-Stratonovich transformation,
$$\Delta\bar{\psi}_{\uparrow}\bar{\psi}_{\downarrow}+\bar{\Delta}\psi_{\uparrow}\psi_{\downarrow},$$
and $\Delta$ plays role of Cooper pairs field. After decoupling, we can safely integrate out fermion fields by introducing Nambu spinor,
$$\Psi=(\psi_{\uparrow}\,\,\bar{\psi}_{\downarrow})^T$$
and its conjugated $\bar{\Psi}$. Integratin over fermion fields gives effective action of theory which is nothing more than GL expansion. Effective action has form
$$\mathrm{tr}\ln(1+\mathcal{G}_0^{-1}\Delta),$$
where
$$\mathcal{G}_0^{-1}=\begin{pmatrix}-\partial_{\tau}+\partial^2/(2m)+\mu & 0\\ 0 & -\partial_{\tau}-\partial^2/(2m)-\mu\end{pmatrix}.$$
Expanding log term and taking trace, you can easily verify that all odd terms have zero contribution.
To see that effective action coincides with GL expansion, it is worth mentioning that we consider only static configurations, so $\Delta$ does not depend on time and effective action will be
$$S_{\text{GL}}=\int d^3x\left[A|\Delta|^2+B|\partial\Delta|^2+C|\Delta|^4\right].$$
Coefficients can be calculated explicitly, gradient term appears if you consider static but spatil inhomogeneous order parameter field.
For me, it is also important that statement about vanishing of odd terms of $\Delta$ in GL expansion does not come from the fact that $n_s\propto |\Delta|^2$. I assume that you derive $n_s$ and I am not sure that quantity $|\Delta|^2$ can be measured.
http://physics.stackexchange.com/questions/243534/invariant-polynomials-of-the-landau-theory-of-phase-transitions-crystal-symmetr
– John M Nov 28 '16 at 19:21