With an increasing number of inverter-interfaced distributed generations (DG) being integrated into distribution systems, significant harmonics could be injected to bring an adverse impact on distribution system operation. These harmonics usually fluctuate due to the variability and intermittency of DG outputs. However, different DGs, depending on their locations and sizes, may present different harmonic contributions in distribution systems, and it is important to identify which DG has the largest influence on a specific bus. In this paper, a Complex Affine arithmetic based Three-phase Harmonic Power Flow method (CATHPF) is proposed to trace the harmonic contributions of individual DGs to distribution systems. The contribution of each DG to harmonic voltage (CDGHV) is calculated based on CATHPF. The CATHPF method is tested on an IEEE 33-bus distribution system. The proposed CATHPF method can quantitatively evaluate harmonic contributions of individual DGs in distribution systems, and in turn guide the effective allocation of harmonic control equipment.