The standard toolkit of operators to probe quanta of geometry in loop quantum gravity consists in area and volume operators as well as holonomy operators. New operators have been defined, in the U(N) framework for intertwiners, which allow to explore the finer structure of quanta of geometry. However these operators do not carry information on the global shape of the intertwiners. Here we introduce dual multipole moments for continuous and discrete surfaces, defined through the normal vector to the surface, taking special care to maintain parametrization invariance. These are raised to multipole operators probing the shape of quantum surfaces. Further focusing on the quadrupole moment, we show that it appears as the Hessian matrix of the large spin Gaussian approximation of coherent intertwiners, which is the standard method for extracting the semi-classical regime of spinfoam transition amplitudes. This offers an improvement on the usual loop quantum gravity techniques, which mostly focus on the volume operator, in the perspective of modeling (quantum) gravitational waves as shape fluctuations waves propagating on spin network states.