Identification of underlying partial differential equations (PDEs) for complex systems remains a formidable challenge. In the present study, a robust PDE identification method is proposed, demonstrating the ability to extract accurate governing equations under noisy conditions without prior knowledge. Specifically, the proposed method combines gene expression programming, one type of evolutionary algorithm capable of generating unseen terms based solely on basic operators and functional terms, with symbolic regression neural networks. These networks are designed to represent explicit functional expressions and optimize them with data gradients. In particular, the specifically designed neural networks can be easily transformed to physical constraints for the training data, embedding the discovered PDEs to further optimize the metadata used for iterative PDE identification. The proposed method has been tested in four canonical PDE cases, validating its effectiveness without preliminary information and confirming its suitability for practical applications across various noise levels.