The aim of this paper is to put the problem of vibroacoustic imaging into the mathematical framework of inverse problems (more precisely, coefficient identification in PDEs) and regularization. We present a model in frequency domain, prove uniqueness of recovery of the spatially varying nonlinearity parameter from measurements of the acoustic pressure at multiple frequencies, and derive Newton as well as gradient based reconstruction methods.