Significance: Quantitative measurement of blood oxygen saturation (sO$_2$) with photoacoustic (PA) imaging is one of the most sought after goals of quantitative PA imaging research due to its wide range of biomedical applications. Aim: A method for accurate and applicable real-time quantification of local sO$_2$ with PA imaging. Approach: We combine multiple illumination (MI) sensing with learned spectral decoloring (LSD); training on Monte Carlo simulations of spectrally colored absorbed energy spectra, in order to apply the trained models to real PA measurements. We validate our combined MI-LSD method on a highly reliable, reproducible and easily scalable phantom model, based on copper and nickel sulfate solutions. Results: With this sulfate model we see a consistently high estimation accuracy using MI-LSD, with median absolute estimation errors of 2.5 to 4.5 percentage points. We further find fewer outliers in MI-LSD estimates compared to LSD. Random forest regressors outperform previously reported neural network approaches. Conclusions: Random forest based MI-LSD is a promising method for accurate quantitative PA oximetry imaging.