This is a straightforward classification problem. The trick is in finding a feature vector for each spectrum. Then we can simply associate each one with the nearest training spectrum in the feature space.
The simplest way to reduce each spectrum to a vector is to represent it as the parameters of a mixture model. Gaussian is the obvious candidate, but a distribution with non-negative support like gamma seems appropriate for your problem too. I'll use Gaussians for illustration since their parameters are easy to estimate by inspection.
We are given the assumption that there are invariably two peaks, so this makes the model very simple: a mean and a variance for each of the two normal densities, plus a relative mixture weight (we can fix one without loss of generality) adding up to five real numbers. In other words, we are mapping each spectrum to a 5D vector.
Your illustration can be roughly approximated as $\mathbf x=[66\; 10 \;103 \;10 \;9]$, where the first four co-ordinates describe the densities and last co-ordinate determines the relative weight of the second density in the mixture. This is what the Gaussian mixture approximation looks like:

I assume you are capable of estimating the parameters of the mixture model with your spectra so I'll skip that part.
Another feature extraction algorithm, as you guessed and cbeleites alluded to below, is PCA. This works by taking the eigendecomposition of the autocovariance of your spectra. The features then are the principal values. The benefit of this approach is that it offers the most compact possible representation of your spectra; i.e., if you allot a fixed number of dimensions to your feature vector, the principal values are the optimal choice to reduce the mean squared error. The drawback is that it works offline (online SVD/PCA variants notwithstanding), while the mixture model algorithm is online.
The final part is calculating the distances. The simplest metric is the Euclidean distance, which equally weights the features. If you want to emphasize certain features over others, such as the mean over the variance, you can use a weighted norm:
$\| x - y \|_p = \left[ \sum_i w_i ( x_i-y_i )^p \right]^{1/p}$, where the weights $w_i$ sum to unity.