Multispectral optoacoustic tomography (MSOT) utilizes multiple wavelengths to illuminate tissue, producing a series of optoacoustic images rich in spectral information. This approach offers a more comprehensive spectral profile compared to conventional optoacoustic techniques. Contrasted with single-wavelength optoacoustic images, the spectral information can be amalgamated with endogenous chromophores or exogenous dyes within biological organisms, thereby unveiling physiological, cellular, and subcellular functions. The development of eigenspectral optoacoustic tomography (eMSOT), grounded in the linear mixture model (LMM), along with its various derivative methods, facilitates label-free imaging of tissue oxygen saturation in deep-seated structures. However, the effectiveness of the LMM may diminish in the presence of multiple scattering effects or inter-substance interactions, thereby impairing the performance of the eMSOT method in heterogeneous tissues. To address this issue, we propose incorporating a nonlinear model to enhance the eMSOT technique, which we refer to as NL-eMSOT (non-linear eMSOT). This model employs the Hadamard product as a nonlinear component of the LMM, effectively characterizing the interactions between photons and both oxygenated and deoxygenated hemoglobin within the near-infrared spectral window. This innovation resolves the nonlinear unmixing problem inherent in optoacoustic imaging. Our approach, validated through numerical simulations, phantom experiments, and in vivo studies, improves the accuracy of quantitative oxygen saturation estimation in heterogeneous tissues by accounting for inter-substance interactions. Consequently, it necessitates the consideration of more complex mixing models to adequately address nonlinear interactions.