Abstract: As crucial nuclear data, the neutron-induced fission yield of ²³⁹Pu plays a vital role in nuclear engineering applications, technological developments, and theoretical studies of fission mechanisms. However, existing experimental data for ²³⁹Pu fission yields across varying neutron energies remain scarce and exhibit significant uncertainties, resulting in insufficient reliable data to support critical studies such as yield-energy relationships, yield-mass distributions, and yield-charge distributions. Therefore, it is necessary to conduct further experimental measurements to obtain higher-precision data. Based on the Xi’an Pulsed Reactor (XAPR) rabbit irradiation facility and quartz-encapsulated plutonium targets, a series of cumulative yield measurements for thermal neutron-induced ²³⁹Pu fission were carried out. A computational formula was derived to calculate the relative yields among various fission products based on their experimentally measured γ-ray peak areas. To account for the differences in half-lives and activities of the target fission products, three experiments with different irradiation durations and powers were designed by optimizing the experimental scheme using a self-developed fission product gamma-ray spectrum simulation program. After irradiation and appropriate cooling, the contents of fission products in the irradiated target were obtained using gamma-ray non-destructive spectrometry, with the measurement duration for a single spectrum ranging from 120 s to 10 d. The gamma-ray detection efficiency curves of non-destructive were calibrated by an equivalent method with several standard planar sources (⁵⁷Co, ¹⁵²Eu, ¹³³Ba, ¹⁸²Ta, and ²⁴Na). To identify the sources and quantify the proportions of interference affecting the target energy peaks during data processing, a self-developed program for fission product gamma full-energy peak interference retrieval was employed. This method was compared with the traditional interference proportion subtraction approach, demonstrating good agreement between the two. However, the former method exhibits a broader scope of applicability. The relative quantitative relationships of the fission products on 41 mass chains such as ⁸⁸Kr, ⁹⁹Mo, ¹³⁸Xe and ¹⁴¹Ce were obtained by processing the measured gamma-ray spectra. Using the ⁹⁹Mo yield value recommended by the ENDF/B-Ⅷ.0 evaluated database as the benchmark, the yields of fission products on other 40 mass chains were determined after some corrections such as influence of non-primary neutron field and non-primary isotope. The uncertainty of the yield measurement results was evaluated, with the uncertainty of the characteristic gamma-ray emission probability being the dominant contributor. Most of the measured values in this work agree with the recommended values in evaluated databases and literature values within the uncertainties. Although the measured values of ¹⁰⁷Rh and ¹²⁷Sb obtained in this work are lower than the literature and evaluated values, the analysis suggests that the present measurement results are more reliable, due to the higher precision of the measurement instrument and the more accurate subtraction method for interference contributions. The sum of the experimentally determined mass chain yields is 149%, and reaches 198% after minor extrapolation and interpolation. To achieve the normalization of mass yields to 200%, it is necessary to conduct further measurements, particularly focusing on short-lived nuclides through the application of rapid chemical separation devices.