You already have a direct and simple solution( flatten first) using basic Mathematica and it works for any depth of nesting.
A different way to do that is to use a Listable function (alias for "==") and flatten the result.
eq[x_, y_] = Equal[x, y]; SetAttributes[eq, Listable];
Flatten[eq[amat, bmat]]
{a11 == b11, a12 == b12, a21 == b21, a22 == b22, a31 == b31,
a32 == b32, a41 == b41, a42 == b42}