When applied to complex engineering system analysis, multi-fidelity Gaussian process regression tends to suffer reduced model accuracy when handling high-dimensional inputs due to the curse of dimensionality. Existing mitigation strategies exhibit limitations such as optimization instability and inadequate feature representation. Targeting this problem, an ensemble deep feature multi-fidelity Gaussian process regression method is proposed. It utilizes an ensemble of deep neural networks to adaptively map high-dimensional inputs to a robust, low-dimensional latent feature space, enhancing representation robustness. A gradient isolation and two-stage training strategy is employed, decoupling the feature extractor pre-training process based on low-fidelity data from the subsequent multi-fidelity Gaussian process regression model construction based on fixed features, circumventing the instability associated with end-to-end optimization in deep fusion models and ensuring robust and efficient training. Finally, the effectiveness of the propsed method is validated through simulations on standard high-dimensional test functions, and its potential for solving practical engineering problems is demonstrated using a case study on equipment range optimization.