Constitutive modeling of heterogeneous hyperelastic materials is still a challenge due to their complex and variable microstructures. We propose a multiscale data-driven approach with a hierarchical learning strategy for the discovery of a generic physics-constrained anisotropic constitutive model for the heterogeneous hyperelastic materials. Based on the sparse multiscale experimental data, the constitutive artificial neural networks for hyperelastic component phases containing composite interfaces are established by the particle swarm optimization algorithm. A microscopic finite element coupled constitutive artificial neural networks solver is introduced to obtain the homogenized stress-stretch relation of heterogeneous materials with different microstructures. And a dense stress-stretch relation dataset is generated by training a neural network through the FE results. Further, a generic invariant representation of strain energy function (SEF) is proposed with a parameter set being implicitly expressed by artificial neural networks (SANN), which describes the hyperelastic properties of heterogeneous materials with different microstructures. A convexity constraint is imposed on the SEF to ensure that the multiscale constitutive model is physically relevant, and the $ℓ_1$ regularization combined with thresholding is introduced to the loss function of SANN to improve the interpretability of this model. Finally, the multiscale model is hierarchically trained, cross-validated and tested using the experimental data of cord-rubber composite materials with different microstructures. The proposed multiscale model provides a convenient and general methodology for constitutive modeling of heterogeneous hyperelastic materials.