Hydrological Simulation Program-Fortran (HSPF) model has many parameters and complex interactions. The traditional parameter optimization is insensitive to the optimization parameters and the optimization algorithm is easy to trap into local problems, which affects the precision and efficiency of optimization. In this paper, a new optimization approach is explored by integrating Qinglong River watershed, parameter sampling, sensitivity analysis, and parameter optimization. Response surface optimization software Design Expert was applied to sample the parameters of 9 HSPF models, and 130 sets of parameter sets were obtained. Multiple quadratic regression models were used to establish the nonlinear relationship between the parameter sets and the efficiency coefficient of nash-sutcliffe (NSE), and the optimal parameters and their dense value ranges were identified by contour lines and response surface. The NSE mean value, maximum value, and minimum value of the response surface optimization parameters as well as the optimized interval reduction rate were all superior to the orthogonal range analysis method. LZSN, INFILT, and AGWRC were extremely sensitive parameters, while DEEPFR was sensitive parameters. The interactions between LZSN and INFILT, INFILT and AGWRC, INFILT and UZSN, and INFILT and IRC had significant impacts on the results. The dense value range of parameters were optimized as follows: LZSN[2.0,2.65], INFILT[0.0,0.475], AGWRC[0.0,0.885], DEEPFR[0.1,0.176], BASETP[0.1,0.120], AGWETP[0,3,0.120], CEPSC[0.6,0.244], UZSN[0.3,1.22], IRC[0.3,0.63]. The response surface method synthesized three aspects, i.e., parameter sampling, parameter sensitivity analysis, and parameter optimization, which considers the nonlinear relationship of parameters, the interaction of parameters, and the optimization accuracy and efficiency, thus opening up a new way for parameter optimization of HSPF model in Qinglong River watershed.