This won't fit in the comment boxes, so let me just add another answer. I think in this case this will be better than reformatting the original question to include the remarks below. 

On second thought I do find Will's model to be a convincing evidence that the logarithmic asymptotic count should likeweise be $d^2 + O(d)$. Below I suggest an explanation for the apparent discrepancy in the two models as far as the distribution of points of very small height is concerned. Taking into account this and the (Schanuel-Schmidt)-Masser-Vaaler theorem, I would think that the statistical parallels between number and function fields are much too precise to think of any answer other than $\log{N(d)} = d^2 + O(d)$. 

Let $n(h,d)$ be the logarithm of the number of points of degree $\leq d$ and logarithmic height $\leq h$. Northcott's bound is $n(h,d) \leq h(d^2+d) + O(d^2)$, which is in particular $O(hd^2)$ when $h$ is bounded away from zero. Its proof is valid also in the function field case, but there, by an extension of Will's argument, we have a precise asymptotic: Assuming as we may that $h \in \frac{\log{q}}{d} \mathbb{Z}^{> 0}$, then in fact $n(h,d) = h(d^2+d) + O(d)$ as $\max(h,d) \to \infty$. In other words: in our function field model, Northcott's bound is always sharp.

For the number field case, Masser and Vaaler, extending previous work of Schanuel and Schmidt, consider the count for $d$ fixed and $h \to \infty$; they prove in particular that $n(h,d) \sim h(d^2+d)$: Northcott's bound is sharp in this situation! (Their result is more precise: it explicitly determines the coefficient of $\exp(h(d^2+d))$, and applies more generally to relative extensions over a fixed number field). 

In $\overline{\mathbb{Q}}^{\times}$, however, in contrast to the situation in our function field model, Northcott's prediction breaks down completely if we take $h = c/d$, where $0 < c < \infty$ is a constant: $\exp(hd^2)$ is then exponential in $d$, whereas we expect $|\{ \alpha \mid h(\alpha) < c/d, \, [\mathbb{Q}(\alpha):\mathbb{Q}] \leq d \}|$ to be polynomially bounded, perhaps even by $\kappa(c)d^2$. (I think we may even be able to prove this; does such a bound appear anywhere in the literature?) 

Here, I believe, is the explanation for this apparent discrepancy (exponential versus polynomial) in the distribution of points of very small height - the torsion points in particular. There are two things to keep in mind. The first is the distinction between general Weil heights and canonical (dynamical, normalized, Neron-Tate) heights: they differ by a bounded amount and therefore they product comparable asymptotics - certainly the same rate of growth - for large heights; whereas the distribution of small points is a very fine intrinsic property of the latter heights. The second point to keep in mind is that our particular function field model here is isotrivial (constant, in fact); this accounts for the profusion of torsion (as well small non-torsion) points. An analogy here is to think of a constant elliptic curve over a complex function field: its points of zero canonical height are all the constant sections (an uncountable set), whereas for non-constant elliptic curves they are just the countable set of torsion points.

In either setting, think of our height $h(\cdot)$ as the canonical dynamical height - the global Tate limit attached by Call and Silverman to a dynamical system on $\mathbb{P}^1$ - of the iteration $z \mapsto z^2$. The difference is that in the function field case, this iteration is isotrivial. Isotriviality has no counterpart in the number field setting, and a more accurate function field model in our problem would be to consider the dynamical height $\hat{h}_f$ attached to a non-isotrivial map of $f : \mathbb{P}^1 \to \mathbb{P}^1$ over $\mathbb{F}_q(t)$; or, if you prefer, the canonical height on a non-isotrivial elliptic curve over $\mathbb{F}_q(t)$. In those cases, as in the arithmetic setting, the number of points of degree $\leq d$ and height $< c/d$ should be similarly bounded polynomially in $d$, and the discrepancy in the count will not occur.

In any case, $\hat{h}_f - h = O_f(1)$, and this implies that for $h \gg_ f 0$ sufficiently large, the logarithmic count $n_f$ for $\hat{h}_f$ still satisfies $n_f(h,d) \asymp h(d^2+d) + O(d)$. (A moment of reflection will show that the proportionality constant should still be $1$, although this could be difficult to prove. Actually even a proof of the existence of a proportionality constant would be interesting, for all we know at this point is that $n_f(h,d)$ is locked between two positive multiples of $h(d^2+d) + O(d)$.) But for any fixed positive $h > 0$, we may construct enough points of height $< h$ by taking inverse images under $f^{\circ N}$ of points of height $< h \cdot (\deg{f})^N$, which will be large enough as soon as $N \gg 0$. Hence, for *any* fixed $h > 0$, we still know that $n_f(h,d) \asymp h(d^2+d) + O(d)$. Surely, the proportionality constant should be $1$, but this might be difficult to prove. 

This now should hold in the number field setting too. 

----------------------------------------------------------------------------------------

Let me finish by adding one more remark about the small point analogy. It is often said that the crux of Lehmer's problem, stating $h(\alpha) > c/d$ unless $\alpha$ is torsion, is in the presence of archimedean places, and that the question is pointlessly trivial in the function field setting. This of course is true for the constant function field model considered in Will's answer, and it is true more generally for polarized dynamical systems with everywhere potential good reduction over a global function field (e.g., abelian varieties with everywhere potential good reduction). However, the real analogy arises when we have degenerations; for instance, for the case of the canonical height on a non-isotrivial elliptic curve over $\mathbb{F}_q(t)$. I believe that Lehmer's problem, in this latter case (a non-isotrivial elliptic), would be just as difficult as the original one for $\mathhbb{G}_m$ over $\bar{\matbb{Q}}$. 

If this is correct then, from a point of view of algebraic dynamics, the essential difficulty in Lehmer's problem is not so much in the failure of the ultrametric triangle inequality in number fields, as it is in the presence of places of degenerating - or chaotic - dynamics (though of course the two points are related). After all, it is a common wisdom in arithmetic geometry that archimedean fibres (or complex spaces) should be regarded as totally degenerate: think of the Tate curve, Mumford's theory of $p$-adic Shottky groups, or the non-triviality of Julia sets in Berkovich analytic spaces (they are always trivial for iterations with good reduction over a non-archimedean local field). 

In view of this, all the subtleties in the distribution of small points in $\mathbb{G}_m(\bar{\mathbb{Q}})$ appear to be present also in a non-isotrivial dynamical system on $\mathbb{P}^1$ over $\mathbb{F}_q(t)$, or in a non-isotrivial elliptic curve over a global function field. But we saw that in such a dynamical system, $n_f(h,d) \asymp h(d^2+d)$ for any fixed $h > 0$ (or for any fixed $d$). So the same should persist in the number field setting too - presumably, with proportionality constant equal to $1$. 

To sum up, as soon as the height is bounded away from zero, Northcott's bound is pretty sharp. This is something I did not expect at the time of asking this question.