Accurate forecasting of electric, cooling and heating loads is an important prerequisite and foundation for the operation scheduling and energy management of integrated energy systems. Leveraging the energy coupling characteristics between multivariate load, this paper constructs a joint prediction model for multivariate load based on bidirectional gated recurrent units (BiGRU) and a progressive layered extraction (PLE) network architecture. Firstly, the meteorological features with high correlation are screened as input features of the model through the maximum information coefficient. Then, the BiGRU network is used to extract the temporal features of the multivariate load time series under the integrated energy system and reconstruct the data in this way. Secondly, for the characteristics of different energy sources that are coupled with each other, the improved progressive hierarchical extraction network structure is proposed, and the coupling features are extracted from the complex and multidimensional data through the multilevel sharing of the feature extraction layer. Finally, by changing the structural parameters of the sub-task tower module, the coupled feature information is differentially fused, and the multiple load prediction results are obtained. The actual example results show that the maximum information coefficient screening method adopted in the article is more suitable for feature selection of meteorological data than the traditional Pearson coefficient screening method, and the proposed BiGRU-PLE multivariate load prediction model can reduce the prediction error by more than 5% compared with the single-task model, and by more than 3% compared with the common multitask model.