In addition to the illuminating answers of Will Sawin and Noam D. Elkies, I want to give a different perspective, coming from analytic number theory over $\mathbb{F}_q[T]$, and ultimaltey relating the question to cubic exponential sums. It does not recover your formula in full, but it does show that $$(\star)\, \sum_{(a_1,a_2,a_3)\in \mathbb{F}_q^3} \binom{N(a_1,a_2,a_3)}{2} = \frac{q^5-22q^4+O(q^{7/2})}{1152}$$ holds when $(q,6)=1$. That is, it recovers the first two terms in your formula. It is also amenable to generalizations. Below, I assume throughout that $(q,6)=1$ even if I do not say so explicitly.
Pointwise results: For any $a_1,a_2,a_3\in \mathbb{F}_q$, let $G$ be the Galois group of the polynomial $f(T)=T^4+a_1T^3+a_2T^2+a_3T+x$ over $\mathbb{F}_q(x)$. S. D. Cohen proved, in 1970, a general result that implies the following: $$N(a_1,a_2,a_3) = \frac{q}{|G|} + O(\sqrt{q})$$ holds as long as the splitting field of $T^4+a_1T^3+a_2T^2+a_3T+x$ does not contain $\mathbb{F}_{q^i}$ for $i>1$. Otherwise, $N(a_1,a_2,a_3)=O(\sqrt{q})$.
Say that $f(T)=T^4+a_1T^3+a_2T^2+a_3T+x$ is 'good' if $G\cong S_4$ and the splitting field does not contain $\mathbb{F}_{q^i}$ ($i>1$). Generically, $f$ is good. This can be made precise: In a later paper, Cohen gave sufficient conditions for a polynomial to be good, see Lemma 1 here. It is also known that 'most' polynomials are Morse polynomials ($f$ is Morse if $f'$ has $\deg f-1$ distinct zeros) and that a Morse polynomial is a good polynomial -- see the discussion and references in section 2.2 of this preprint of Kurlberg and Rosenzweig. These results can be used to show $$\sum_{(a_1,a_2,a_3)\in \mathbb{F}_q^3} \binom{N(a_1,a_2,a_3)}{2} \sim q^3 \binom{\frac{q}{24}}{2} \sim \frac{q^5}{1152}.$$
Character approach: Let $\mathcal{M}_q$ be the set of monic polynomials in $\mathbb{F}_q[T]$. Say that $f_1,f_2 \in \mathcal{M}_q$ are '$j$-equivalent' if they have the same first $j$ next-to-leading coefficient (the $j$-th next to leading coefficient of $f$ is the coefficient of $T^{\deg f- j}$ in $f$; if $j>\deg f$ this is $0$). One can detect the condition that $f_1$ and $f_2$ are $j$-equivalent using a variant of Dirichlet characters: let $G(j)$ be the group of characters $$\chi' \colon ((1+T^{-1}\mathbb{F}_q[T^{-1}])/T^{-j-1}\mathbb{F}_q[T^{-1}])^{\times} \to \mathbb{C}^{\times}.$$ There are $q^j$ such characters, and any such character can be extended to a multiplicative function $\chi\colon \mathcal{M}_q \to \mathbb{C}^{\times}$ by $\chi(f) = \chi'(T^{-\deg f} f(T) \bmod T^{-j-1})$. It is easy to check that $$q^{-j}\sum_{\chi \in G(j)} \chi(f_1)\overline{\chi(f_2)} = \mathbf{1}_{f_1 \text{ is equivalent to }f_2}.$$ (If this looks funny, a different construction of such $\chi$-s proceeds taking Dirichlet characters $\chi'$ modulo $T^{j+1}$, which are even (constant on scalars) and then letting $\chi(f):=\chi'( f(1/T)T^{\deg f}/f(0))$ and $\chi(T)=1$.)
This is incredibly relevant because $N(a_1,a_2,a_3)$ is a sum over polynomials of degree $4$ that are $3$-equivalent to $T^4+a_1T^3+a_2T^2+a_1T$ (satisfying some factorization conditions). Letting $$S = \{ f \in \mathcal{M}_q: \deg f = 4, \, f \text{ is a product of 4 linear factors over } \mathbb{F}_q \}$$ we can write $$N(a_1,a_2,a_3) = q^{-3} \sum_{\chi \in G(3)}\overline{\chi( T^4+a_1T^3+a_2T^2+a_3T) }\sum_{f \in S} \chi(f) .$$ By orthogonality relations of the characters in $G(3)$, one obtains $$\sum_{(a_1,a_2,a_3) \in \mathbb{F}_q^3} N^2(a_1,a_2,a_3) = q^{-3} \sum_{\chi \in G(3)}\left| \sum_{f \in S} \chi(f)\right|^2.$$ For the trivial character $\chi_0 \in G(3)$, $\sum_{f \in S} \chi_0(f) =|S|=\binom{q}{4}$. Moreover, $\sum_{(a_1,a_2,a_3) \in \mathbb{F}_q^3} N(a_1,a_2,a_3) = |S| = \binom{q}{4}$. It follows that $$\sum_{(a_1,a_2,a_3)\in \mathbb{F}_q^3} \binom{N(a_1,a_2,a_3)}{2} = \frac{1}{2} \left[q^{-3} \sum_{\chi_0\neq \chi \in G(3)}\left| \sum_{f \in S} \chi(f)\right|^2 + q^{-3}\binom{q}{4}^2 - \binom{q}{4}\right].$$ Using Newton's identities, one can write $$\sum_{f \in S} \chi(f) = \frac{1}{4!} \left( p_1^4 - 6p_1^2 p_2+3 p_2^2+8p_1p_3 - 6 p_4\right)$$ where $$p_k := \sum_{a \in \mathbb{F}_q} \chi^k(T+a).$$ By Weil's RH, each $p_k$ is $O(\sqrt{q})$ when $\chi\neq\chi_0$, which gives $\sum_{f \in S} \chi(f) = O(q^2)$, implying $$\sum_{(a_1,a_2,a_3)\in \mathbb{F}_q^3} \binom{N(a_1,a_2,a_3)}{2} =\frac{1}{2}q^{-3}\binom{q}{4}^2 + O(q^4) =\frac{q^5+O(q^4)}{1152}.$$ Relation to random matrix theory: Let us refine the above result. The total contribution of $\chi_0\neq \chi \in G(2)$ is $O(q^3)$ and we focus on $\chi \in G(3)\setminus G(2)$. By Newton's identities as above, analyzing $\sum_{f \in S} \chi(f)$ is the same as analyzing $(\sum_{a \in \mathbb{F}_q} \chi(T+a))^4$ (up to negligible error), which can be formally written as $(-q^{\frac{1}{2}} \mathrm{Tr}(\Theta_{\chi}))^4$ for $\Theta_{\chi} \in U(2)$ if $\chi$ is primitive (recall $\sum_{a \in \mathbb{F}_q} \chi(T+a)$ is the linear term in an L-function $L(u,\chi)$). Using a result of Katz, building on Deligne's equidistribution theorem, it follows that $$q^{-3} \sum_{\chi \in G(3)\setminus G(2)}\left| q^{-\frac{1}{2}}\sum_{a \in \mathbb{F}_q} \chi(T+a)\right|^2 = \int_{U(2)} |\mathrm{Tr}(U)|^4{\rm d}U + O(q^{-1/2})=14+O(q^{-1/2})$$ and $(\star)$ follows after some algebra. The random matrix integral was evaluated using Theorem 1.1 of Rains.
Relation to cubic exponential sums: The characters in $G(R_3)$ can be parametrized explicitly -- it is $\cong \mathbb{F}_q^3$). Indeed, if we fix a nontrivial character $\psi \colon \mathbb{F}_q \to \mathbb{C}^*$ then for any $x_1,x_2,x_3\in \mathbb{F}_q$ we can define $$\chi_{x_1,x_2,x_3}'(T^n+a_1T^{n-1}+a_2 T^{n-2}+a_3T^{n-3}+\ldots) = \psi( \sum_{i=1}^{3}x_i h_i(a_1,a_2,a_3))$$ where $$h_1(a_1,a_2,a_3)=a_1,\, h_2(a_1,a_2,a_3) = a_2 - \frac{a_1^2}{2},\, h_3(a_1,a_2,a_3) = a_3-a_1a_2 + \frac{a_1^3}{3}.$$ In particular, for the $\chi$ that arises from $\chi'_{x_1,x_2,x_3}$ we have $$\sum_{a \in \mathbb{F}_q} \chi(T+a) = \sum_{a \in \mathbb{F}_q} \psi\left( x_1 a - x_2 \frac{a^2}{2} + x_3 \frac{a^3}{3}\right),$$ which is a cubic exponential sums. High moments of such sums were computed by Livné. In a separate paper he also related them to Galois representations, which appeared in Will Sawin's answer.