Estimation of the bulk optical properties of turbid samples from spatially resolved reflectance measurements remains challenging, as the relation between the bulk optical properties and the acquired spatially resolved reflectance profiles is influenced by wavelength-dependent properties of the measurement system. The resulting measurement noise is apparent in the estimation of the bulk optical properties. In this study, a constrained inverse metamodeling approach is proposed to overcome these problems. First, a metamodel has been trained on a set of intralipid phantoms covering a wide range of optical properties to link the acquired spatially resolved reflectance profiles to the respective combinations of bulk optical properties (absorption coefficient and reduced scattering coefficient). In this metamodel, the wavelength (500 - 1700 nm) is considered as a third input parameter for the model to account for the wavelength dependent effects introduced by the measurement system. Secondly, a smoothness constraint on the reduced scattering coefficient spectra was implemented in the iterative inverse estimation procedure to robustify it against measurement noise and increase the reliability of the obtained bulk absorption and reduced scattering coefficient spectra. As the estimated values in some regions may be more reliable than others, the difference between simulated and measured values as a function of the evaluated absorption and scattering coefficients was combined in a 2D cost function. This cost function was used as a weight in the fitting procedure to find the parameters of the µ(s)' function giving the lowest cost over all the wavelengths together. In accordance with previous research, an exponential function was considered to represent the µ(s)' spectra of intralipid phantoms. The fitting procedure also provides an absorption coefficient spectrum which is in accordance with the measurements and the estimated parameters of the exponential function. This robust inverse estimation algorithm was validated on an independent set of intralipid® phantoms and its performance was also compared to that of a classical single-wavelength inverse estimation algorithm. While its performance in estimating µ(a) was comparable (R2 of 0.844 vs. 0.862), it resulted in a large improvement in the estimation of µ(s)' (R2 of 0.987 vs. 0.681). The change in performance is more apparent in the improvement of RMSE of µ(s)', which decreases from 10.36 cm(-1) to 2.10 cm(-1). The SRS profiles change more sensitively as a function of µ(a). As a result, there is a large range of µ(s)' and a small range of µa resulting in a good fit between measurement and simulation. The robust inverse estimator incorporates information over the different wavelengths, to increase the accuracy of µ(s)'estimations and robustify the estimation process.