This is a nice question. In my view, the choice of test function is motivated by the idea of taking some positive quantity that is a supersolution to an elliptic equation, and looking at the equation for its logarithm. The new equation contains a useful term that is quadratic in the gradient. This idea is pervasive in geometry and elliptic PDEs, and some examples are below.
(1) The basic case to consider is that $u$ is positive and superharmonic. Then $v := -\log u$ satisfies $|\nabla v|^2 \leq \Delta v$, which gives local bounds on the integral of $|\nabla v|^2$ (independent of $v$) after multiplying by standard cutoffs and integrating by parts. This is already enough to prove the Harnack inequality for harmonic functions in two dimensions, since the Dirichlet energy controls oscillation for functions that satisfy the maximum and minimum principle.
(2) In your context, the choice of test function $H(u)$ satisfies $$a^{ij}\partial_iu\partial_j(H(u)) = a^{ij}\partial_i(V(u))\partial_j(V(u)),$$ where $V(u) = c_1\log(c_2 - u)$ with $c_2 - u$ positive. I view the estimate as coming from integrating the equation for $V(u)$.
To illustrate how this works in a simple context, assume that $u \in C^2_0(B_1)$ satisfies $\Delta u \geq -A$, and that $u \leq M$. Then for $w := M+A-u > 0$ we have that $v := \log(M+A)-\log(w)$ is compactly supported and satisfies $|\nabla v|^2 \leq 1 + \Delta v$. Thus the integral of $|\nabla v|^2$ (hence $v^{2^*}$) is bounded in terms of the volume of the domain.
(3) The Bombieri-De Giorgi-Miranda interior gradient estimate for a solution $u$ to the minimal surface equation is based on the fact that the vertical component $\nu^{n+1}$ of the unit normal to the graph of $u$ is positive and superharmonic (on the graph). The proof uses the equation for $v := -\log(\nu^{n+1})$, which just as above contains a useful term quadratic in $|\nabla v|$.
(4) The Li-Yau proof of the Harnack inequality for a harmonic function $u$ is obtained by looking at the quantity $w := |\nabla (-\log u)|^2$. The key is that $w$ solves a differential inequality with the powerful term $\frac{2}{n}w^2$, which allows one to bound $w$ from above locally by a universal constant independent of $w$.
I am sure there are many other interesting examples, and I am not sure where the first instances of the "log trick" appeared. One final remark is that the estimate (4.9) can also quickly be inferred using the properties of the Green's function $G$ for uniformly elliptic operators (namely, $G \in L^p$ for $p < \frac{n}{n-2}$ and $\nabla G \in L^p$ for $p < \frac{n}{n-1}$, just like the Laplace case).