Using $e^+e^-$ annihilation data corresponding to an integrated luminosity of 2.93 $\rm fb^{-1}$ taken at the center-of-mass energy $\sqrt{s}=3.773$~GeV with the BESIII detector, a joint amplitude analysis is performed on the decays $D^0\to\pi^+\pi^-\pi^+\pi^-$ and $D^0\to\pi^+\pi^-\pi^0\pi^0$(non-$\eta$). The fit fractions of individual components are obtained, and large interferences among the dominant components of $D^{0}\to a_{1}(1260)\pi$, $D^{0}\to\pi(1300)\pi$, $D^{0}\to\rho(770)\rho(770)$ and $D^{0}\to2(\pi\pi)_{S}$ are found in both channels. With the obtained amplitude model, the $CP$-even fractions of $D^0\to \pi^+\pi^-\pi^+\pi^-$ and $D^0\to\pi^+\pi^-\pi^0\pi^0$(non-$\eta$) are determined to be $(75.2\pm1.1_{\rm stat.}\pm1.5_{\rm syst.})\%$ and $(68.9\pm1.5_{\rm stat.}\pm 2.4_{\rm syst.})\%$, respectively. The branching fractions of $D^0\to \pi^+\pi^-\pi^+\pi^-$ and $D^0\to\pi^+\pi^-\pi^0\pi^0$(non-$\eta$) are measured to be $(0.688\pm0.010_{\rm stat.}\pm 0.010_{\rm syst.})\%$ and $(0.951\pm0.025_{\rm stat.}\pm 0.021_{\rm syst.})\%$, respectively. The amplitude analysis provides an important model for binning strategy in the measurements of the strong phase parameters of $D^0 \to 4\pi$ when used to determine the CKM angle $\gamma (\phi_{3})$ via the $B^{-}\to D K^{-}$ decay.