Despite advances in seismology and computing, the ability to image subsurface volcanic environments is poor, limiting our understanding of the overall workings of volcanic systems. This is related to substantive structural heterogeneities which strongly scatters seismic waves obscuring the ballistic arrivals normally used in seismology for wave velocity determination. Here we address this constraint by, using a deep learning approach, a Fourier neural operator (FNO), to model and invert seismic signals in volcanic settings. The FNO is trained using 40,000+ simulations of elastic wave propagation through complex volcano models, and includes the full scattered wavefield. Once trained, the forward network is used to predict elastic wave propagation and is shown to accurately reproduce the seismic wavefield. The FNO is also trained to predict heterogeneous velocity models given a limited set of input seismograms. It is shown to capture details of the complex velocity structure that lie far outside the ability of current methods available in volcano imagery.