This paper addresses the problem of tracking multiple extended targets in three-dimensional space. We propose the Gaussian process Gaussian mixture probability hypothesis density (GP-PHD) filter, which is capable of tracking multiple extended targets with complex shapes in the presence of clutter. Our approach combines the Gaussian process regression measurement model with the probability hypothesis density filter to estimate both the kinematic state and the shape of the targets. The shape of the extended target is described by a 3D radial function and is estimated recursively using the Gaussian process regression model. Furthermore, we transform the recursive Gaussian process regression problem into a state estimation problem by deriving a state space model such that the estimation of the extent can be integrated into the kinematic part. We derive the predicted likelihood function of the PHD filter and provide a closed-form Gaussian mixture implementation. To evaluate the performance of the proposed filter, we simulate a typical extended target tracking scenario and compare the GP-PHD filter with the traditional Gamma Gaussian Inverse-Wishart PHD (GGIW-PHD) filter. Our results demonstrate that the proposed algorithm outperforms the GGIW-PHD filter in terms of estimating both kinematic states and shape. We also investigate the impact of the measurement rates on both filters; it is observed that the proposed filter exhibits robustness across various measurement rates, while the GGIW-PHD filter suffers under low-measurement-rate conditions.