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 ultimately showing the question is equivalent to one about moments of **cubic exponential sums**. The computations below won't recover your formula in full, but they do 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, they recover the first two terms in your formula. With more work this perspective should recover the full formula. The advantage of this perspective is that it is amenable to generalizations. From now on I assume 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][1] 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][2]. 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 work of Kurlberg and Rosenzweig][3]. 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 coefficients (the $i$-th next to leading coefficient of $f$ is the coefficient of $T^{\deg f- i}$ in $f$; if $i>\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 (abusing notation) 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})$. We abuse notation again and use $G(j)$ for the set of these functions on $\mathcal{M}_q$. 

(If this looks funny, a different construction of such $\chi$-s  proceeds by taking the even Dirichlet characters $\chi'$ modulo $T^{j+1}$ and then letting $\chi(f):=\chi'( f(1/T)T^{\deg f}/f(0))$ and $\chi(T)=1$. A 3rd construction is given at the end of this answer.)

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$ (and satisfying some factorization conditions). To mage use of $G(j)$ we employ two orthogonality relations:
$$(O1)\, q^{-j}\sum_{\chi \in G(j)} \chi(f_1)\overline{\chi(f_2)} = \mathbf{1}_{f_1 \text{ is }j\text{-equivalent to }f_2}$$
for any $f_1,f_2 \in \mathcal{M}_q$ and 
$$(O2)\, q^{-j}\sum_{f \in X_j} \chi_1(f)\overline{\chi_2(f)} = \mathbf{1}_{\chi_1=\chi_2} $$
where $\chi_1,\chi_2 \in G(j)$ and $X_j$ is a set of $q^j$ polynomials, none of which is $j$-equivalent to the other. Letting
$$S = \{ f \in \mathcal{M}_q: \deg f = 4, \, f \text{ is a product of distinct 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 $(O1)$ and $(O2)$, 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 Riemann Hypothesis, 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 written as $(-q^{\frac{1}{2}} \mathrm{Tr}(\Theta_{\chi}))^4$ for some $\Theta_{\chi} \in U(2)$ ($\sum_{a \in \mathbb{F}_q} \chi(T+a)$ is just the linear term in an $L$-function $L(u,\chi)$ which has two zeros on $|u|=q^{-1/2}$). Using [Theorem 1.2 of Katz][4], 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))^4\right|^2 = \int_{U(2)} |\mathrm{Tr}(U)|^8{\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][5].

**Relation to cubic exponential sums:** The characters in $G(3)$ can be parametrized explicitly and $G(3)\cong \mathbb{F}_q^3$. Indeed, if we fix a nontrivial character $\psi \colon \mathbb{F}_q \to \mathbb{C}^*$ we can define an isomorphism $\mathbb{F}_q^3 \to G(3)$ by sending $x_1,x_2,x_3\in \mathbb{F}_q$ to
$$\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, 
$$\sum_{a \in \mathbb{F}_q} \chi_{x_1,x_2,x_3}(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 sum. Birch computed the first 4 even moments of cubic sums in 1968 (see [last page of his paper][6]), obtaining
$$\sum_{a \in \mathbb{F}_p^{\times},\, b \in \mathbb{F}_p}\left| p^{-1/2}\sum_{x \in \mathbb{F}_p}e_p(ax^3+bx)\right|^{2R} = \frac{(2R-1)!}{(R-1)!(R+1)!}(p-1) (2p-R+1)$$
for $p>3$ and $1\le R \le 4$ where $e_p(x) =e^{2\pi i x/p}$. Higher moments were [computed by Livné][7]. For $R=4$ this is $\sim 14p^2$, consistent with Katz's work. 

In a separate paper Livné [related cubic sums to Galois representations][8], which appeared in Will Sawin's answer.


  [1]: https://zbmath.org/0209.36001
  [2]: https://zbmath.org/0253.12019
  [3]: https://arxiv.org/abs/1802.01215
  [4]: https://zbmath.org/1328.13028
  [5]: https://zbmath.org/0885.05112
  [6]: https://zbmath.org/0183.25503
  [7]: https://zbmath.org/0602.10027
  [8]: https://zbmath.org/0621.14019