Nonlinear optimization methods are very effective in finding models that predict optimal waveform fits to observed seismograms, provided the number of model parameters is small. To reduce the size of the large-scale two- or three-dimensional seismic modelling problem, path integrals over seismic velocities may be defined as model parameters instead of the seismic velocities themselves. This partitions the large-scale optimization problem into a number of independent small-scale problems and reduces the computational effort by several orders of magnitude. The path integral parameters are transformed such that their correlation matrix becomes diagonal. A linear inversion performs the mapping from the space of path integral parameters into the model space of seismic velocities and makes Backus-Gilbert resolution analysis possible. The method has been applied to body wave and surface wave data from the broadband portable Network of Autonomously Recording Seismographs (NARS) in western Europe. Using a surface wave ray approximation to generate synthetics by mode summation, waveforms were fitted for frequencies as high as 50 mHz. Variations in S velocity as large as 6¿1% have been established with a vertical resolution of 100--150 km and a horizontal resolution as small as 400 km in parts of the model. Although a low velocity layer is present under the southern part of the NARS array, such layer is very weak or absent north of the Paris Basin. Horizontal velocity gradients are well-resolved down to 400 km depth. ¿ American Geophysical Union 1990