Harmonics and three-phase imbalance are two common disturbances in power systems that interact with and affect each other, leading to more complex impacts on power system performance. To comprehensively and accurately analyze the distribution and interaction mechanisms of harmonics and three-phase imbalance, this study investigated the joint emission level of harmonic and three-phase unbalanced disturbance sources (H-TU-DS). This paper proposes a joint probabilistic emission level modeling approach. The power, harmonic, and three-phase unbalanced data of H-TU-DS are multidimensional, and the feature dataset is filtered by calculating multiple correlation coefficients to achieve dimensionality reduction without losing data information. Considering the different operating characteristics of the H-TU-DS, an improved fuzzy C-means (FCM) algorithm is applied to classify operating conditions, enabling a comprehensive and accurate analysis of the joint emission level. Based on the nonparametric kernel density function, joint probability modeling of H-TU-DS is performed using an adaptive bandwidth improvement method, effectively enhancing the accuracy of the modeling process. Finally, the effectiveness and accuracy of the proposed methods for feature dataset filtering, operating condition classification, and joint probability modeling are demonstrated through comparisons with both simulated and actual data. The proposed model facilitates the prediction and assessment of the impact of H-TU-DS on power quality after grid connection, enabling the effective prevention of power-quality problems through advanced treatment.