The traditional interpolation might cause smoothing effect on the geochemical elements’ concentration due to its moving weighted average properties. Meanwhile, it cannot describe the elements’ spatial variability if the data’s distribution is sparse, which could cause uncertainties when inferring the information at the unsampled locations. Since the Multiple-Point Geostatistics (MPS) is a stochastic simulation based on the pattern’s statistical characteristics which conducted by spatial multiple points, it can reduce the smoothing effect and quantify the spatial variabilities and uncertainties. Meanwhile, the Local Singularity Analysis (LSA) could identify the weak geochemical anomalies more effectively due to its statistical self-similarity. Therefore, we proposed a new method combing the advantages of both MPS and LSA for improving the identification of the weak geochemical anomalies and quantifying the uncertainty distribution. We firstly treated the geochemical raster data as the Training Images (TIs) and then the Direct Sampling (DS) algorithm was used to perform elements spatial distribution. The LSA and information fusion were finally adopted to obtain the probability map of geochemical anomalies. We verified our methods by utilizing the 1:200,000 stream sediments geochemical data, and the results show that these identified anomalous fields could overcame the smoothing effect and quantify the distribution uncertainties on different simulation results. According to the ROC curves’ analysis, this method has a better performance than the Kriging-based ones and its anomalous fields were more consistent with the determined Fe deposits, which can help us have a clearer understanding for the metallogenic backgrounds and the element distribution patterns.