Phase diversity phase retrieval (PDPR) has been a popular technique for quantitatively measuring wavefront errors of optical imaging systems by extracting the phase information from several designated intensity measurements. As the problem is inverse and non-convex in general, the accuracy and robustness of most such algorithms rely greatly on the initial conditions. In this work, we propose a new strategy that combines Limited-Memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) with the initial points generated by k-means clustering method and three various channels to improve the overall performance. Experimental results show that, for 500 different phase aberrations with root mean square (RMS) value bounded within [0.2λ, 0.3λ], the minimum, the maximum and the mean RMS residual errors reach 0.017λ, 0.066λ and 0.039λ, respectively, and 84.8% of the RMS residual errors are less than 0.05λ. We have further investigated and analyzed the proposed method in details to quantitatively demonstrate its performance: statistical results reveal that our proposed PDPR with k-means clustering enhanced method has excellent robustness in terms of initial points and other influential factors, and the accuracy can outperform its counterpart methods such as classic L-BFGS and modified BFGS.