The Turing patterns of reaction-diffusion equations defined over a square region are more complex because of the D4-symmetry of the spatial region. This leads to the occurrence of multiple equivariant Turing bifurcations. In this paper, taking the FHN model as an example, we give a explicit calculation formula of normal form for the simple and double Turing bifurcation of the reaction-diffusion equation with Dirichlet boundary conditions and defined on a square space, and we also obtain a method for the calculation of the existence of spatially inhomogeneous steady-state solutions. This paper provides a theoretical basis for exploring and predicting the pattern formation of spatial multimode interaction.