Aiming at the problem of dynamic model reduction of mistuned blisks, this study proposes a model reduction method based on the proper orthogonal decomposition (POD) in the frequency domain and establishes the reduced-order model (ROM). The mistuning amount of blisk is characterized as generalized mistuned load, and then a mistuned load interpolation algorithm is introduced to solve the forced vibration response. POD ROM is used to analyze the modal characteristics and vibration response characteristics of a mistuned blisk with 20 sectors, and the results are compared with the full-order model (FOM) and the subset nominal mode (SNM) reduced model. The results show that the model reduction method proposed in this paper can improve the calculation efficiency while ensuring calculation accuracy. In the aspect of modal characteristics of the mistuned blisk, its calculation accuracy is better than the SNM model reduction method, the calculation error of natural frequency is less than 0.05%, and the spatial distribution characteristics of mistuned modes are closer to the results of FOM. In terms of response characteristics, the forced response accuracy is affected by the integration time step and relaxation factor, and the steady-state amplitude error in the time domain can be controlled below 5%.