Retrieving　multi-temporal　and　large-scale　thermohaline　structure　information　of　the　interior　of　the　global　ocean　based　on　surface　satellite　observations　is　important　for　understanding　the　complex　and　multidimensional　dynamic　processes　within　the　ocean.　This　study　proposes　a　new　ensemble　learning　algorithm,　extreme　gradient　boosting　(XGBoost),　for　retrieving　subsurface　thermohaline　anomalies,　including　the　subsurface　temperature　anomaly　(STA)　and　the　subsurface　salinity　anomaly　(SSA),　in　the　upper　2000　m　of　the　global　ocean.　The　model　combines　surface　satellite　observations　and　in　situ　Argo　data　for　estimation,　and　uses　root-mean-square　error　(RMSE),　normalized　root-mean-square　error　(NRMSE),　and　R-2　as　accuracy　evaluations.　The　results　show　that　the　proposed　XGBoost　model　can　easily　retrieve　subsurface　thermohaline　anomalies　and　outperforms　the　gradient　boosting　decision　tree　(GBDT)　model.　The　XGBoost　model　had　good　performance　with　average　R-2　values　of　0.69　and　0.54,　and　average　NRMSE　values　of　0.035　and　0.042,　for　STA　and　SSA　estimations,　respectively.　The　thermohaline　anomaly　patterns　presented　obvious　seasonal　variation　signals　in　the　upper　layers　(the　upper　500　m);　however,　these　signals　became　weaker　as　the　depth　increased.　The　model　performance　fluctuated,　with　the　best　performance　in　October　(autumn)　for　both　STA　and　SSA,　and　the　lowest　accuracy　occurred　in　January　(winter)　for　STA　and　April　(spring)　for　SSA.　The　STA　estimation　error　mainly　occurred　in　the　El　Nino-Southern　Oscillation　(ENSO)　region　in　the　upper　ocean　and　the　boundary　of　the　ocean　basins　in　the　deeper　ocean;　meanwhile,　the　SSA　estimation　error　presented　a　relatively　even　distribution.　The　wind　speed　anomalies,　including　the　u　and　v　components,　contributed　more　to　the　XGBoost　model　for　both　STA　and　SSA　estimations　than　the　other　surface　parameters;　however,　its　importance　at　deeper　layers　decreased　and　the　contributions　of　the　other　parameters　increased.　This　study　provides　an　effective　remote　sensing　technique　for　subsurface　thermohaline　estimations　and　further　promotes　long-term　remote　sensing　reconstructions　of　internal　ocean　parameters.