  Dust emission is the main foreground to Cosmic Microwave Background (CMB)
polarization. Its statistical characterization must be derived from the
analysis of observational data because the precision required for a reliable
component separation is far greater than currently achievable with physical
models of the turbulent magnetized interstellar medium. This letter takes a
significant step towards this goal by proposing a method that retrieves
non-Gaussian statistical characteristics of dust emission from noisy Planck
polarization observations at 353 GHz. We devise a statistical denoising method
based on the Wavelet Phase Harmonics (WPH) statistics, which characterize the
coherent structures in non-Gaussian random fields and define a generative model
of the data. The method is validated on mock data combining a dust map from a
magnetohydrodynamic simulation and Planck noise maps. The denoised map
reproduces the true power spectrum down to scales where the noise power is an
order of magnitude larger than that of the signal. It remains highly correlated
to the true emission and retrieves some of its non-Gaussian properties. Applied
to Planck data, the method provides a new approach towards building a
generative model of dust polarization that will characterize the full
complexity of the dust emission.
