The increasing penetration of wind power brings great uncertainties into power systems, which poses challenges to system planning and operation. This paper proposes a novel probabilistic load flow (PLF) method based on clustering technique to handle large fluctuations from large-scale wind power integration. The traditional cumulant method (CM) for PLF is based on the linearization of load flow equations around the operating point, therefore resulting in significant errors when input random variables have large fluctuations. In the proposed method, the samples of wind power and loads are first generated by the inverse Nataf transformation and then clustered using an improved K-means algorithm to obtain input variable samples with small variances in each cluster. With such pre-processing, the cumulant method can be applied within each cluster to calculate cumulants of output random variables with improved accuracy. The results obtained in each cluster are combined according to the law of total probability to calculate the final cumulants of output random variables for the whole samples. The proposed method is validated on modified IEEE 9-bus and 118-bus test systems with additional wind farms. Compared with the traditional CM, 2m+1 point estimate method (PEM), Monte Carlo simulation (MCS) and Latin hypercube sampling (LHS) based MCS, the proposed method can achieve a better performance with consideration of both computational efficiency and accuracy.