This work is concerned with discovering the governing partial differential equation (PDE) of a physical system. Existing methods have demonstrated the PDE identification from finite observations but failed to maintain satisfying results against noisy data, partly owing to suboptimal estimated derivatives and found PDE coefficients. We address the issues by introducing a noise-aware physics-informed machine learning (nPIML) framework to discover the governing PDE from data following arbitrary distributions. We propose training a couple of neural networks, namely solver and preselector, in a multi-task learning paradigm, which yields important scores of basis candidates that constitute the hidden physical constraint. After they are jointly trained, the solver network estimates potential candidates, e.g., partial derivatives, for the sparse regression to initially unveil the most likely parsimonious PDE, decided according to information criterion. Denoising physics-informed neural networks (dPINNs), based on Discrete Fourier Transform (DFT), is proposed to deliver the optimal PDE coefficients respecting the noise-reduced variables. Extensive experiments on five canonical PDEs affirm that the proposed framework presents a robust and interpretable approach for PDE discovery, leading to a new automatic PDE selection algorithm established on minimization of the information criterion decay rate.