Ocean sound speed is important for underwater acoustic applications, such as communications, navigation and localization, where the assumption of uniformly distributed sound speed profiles (SSPs) is generally used and greatly degrades the performance of underwater acoustic systems. The acquisition of SSPs is necessary for the corrections of the sound ray propagation paths. However, the inversion of SSPs is challenging due to the intricate relations of interrelated physical ocean elements and suffers from the high costs of calculations and hardware deployments. This paper proposes a novel sound speed inversion method based on multi-source ocean remote sensing observations and machine learning, which adapts to large-scale sea regions. Firstly, the datasets of SSPs are generated utilizing the Argo thermohaline profiles and the empirical formulas of the sound speed. Then, the SSPs are analyzed utilizing the empirical orthogonal functions (EOFs) to reduce the dimensions of the feature space as well as the computational load. Considering the nonlinear regression relations of SSPs and the observed datasets, a general framework for sound speed inversion is formulated, which combines the designed machine learning models with the reduced-dimensional feature representations, multi-source ocean remote sensing observations and water temperature data. After being well trained, the proposed machine learning models realize the accurate inversion of the targeted ocean region by inputting the real-time ocean environmental data. The experiments verify the advantages of the proposed method in terms of the accuracy and effectiveness compared with conventional methods.