Assume b1b2 ≠ 0 and solve the linear system for (a1², a2²). Add the two equations, so you get a1² + a2² = a3² b3 / b1, and conclude that b1 divides a3 (using that the b's are pairwise coprime, squarefree). Similarly subtract the two equations to get b2 divides a4. Use this to conclude that b1 divides a2 as well (e.g. by substituting into the original equations, or working in Z/b1 Z), and b2 divides a1 as well. Proceed to conclude b2 divides a3, etc. until you have that all the a's are multiples of b1 b2.

Now looking back at the solution of the linear system for (a1², a2²), we still have a pesky factor 1/2 (coming from the determinant), but at least we see that a two-parameter family of solutions to the original problem is given by a3 = b1² b2 b4 x, a4 = b1 b2² b3 y with x,y of the same parity, and a1, a2 given by the solution to the linear system; explicitly a1² = (1/2) b1² b2² b3 b4 (x + y), a2² = (1/2) b1² b2² b3 b4 (x - y).