The high penetration of photovoltaic (PV) generation complicates the transient characteristics of distribution networks, leading to insufficient feature extraction and weak topological correlation in traditional high-impedance fault (HIF) diagnosis and location methods, thereby degrading accuracy. To address these challenges, a novel HIF diagnosis and location method for PV-integrated distribution networks is proposed based on multi-feature fusion via zebra optimization algorithm (ZOA)-optimized multivariate variational mode decomposition (MVMD), combined with Teager-Kaiser energy operator (TKEO), and a spatiotemporal graph convolutional network (STGCN). First, MVMD processes multivariate signals to effectively fuse multi-dimensional data and fully extract fault features. Meanwhile, ZOA is used to optimize the MVMD parameters, further enhancing feature extraction performance. Second, TKEO is utilized to enhance the highest-frequency intrinsic mode function decomposed by MVMD, enabling the capture of transient energy variations. Finally, a multi-feature fusion vector set is constructed and fed into the STGCN, which employs long short-term memory layers to extract temporal dynamic features and integrates graph convolutional networks to explore spatial topological relationships among nodes, achieving joint spatiotemporal feature modeling. Simulation tests conducted on the IEEE 33-node system demonstrate that, compared to traditional methods, the proposed approach achieves higher accuracy in HIF diagnosis and location for PV distribution networks.