Shale gas has been playing an increasingly important role in meeting global energy demands. The heterogeneity of the pore structure in organic-rich shales greatly affects the adsorption, desorption, diffusion and flow of gas. The pore size distribution (PSD) is a key parameter of the heterogeneity of the shale pore structure. In this study, the Neimark-Kiselev (N-K) fractal approach was applied to investigate the heterogeneity in the PSD of the lower Silurian organic-rich shales in South China using low-pressure N₂ adsorption, total organic carbon (TOC) content, maturity analysis, X-ray diffraction (XRD) and field emission scanning electron microscopy (FE-SEM) measurements. The results show that (1) the fractal dimension DN-K obtained by N-K theory better represents the heterogeneity of the PSD in shale at an approximately 1-100 nm scale. The DN-K values range from 2.3801 to 2.9915, with a mean of 2.753. The stronger the PSD heterogeneity is, the higher the DN-K value in shale is. (2) The clay-rich samples display multimodal patterns at pore sizes greater than 20 nm, which strongly effect the PSD heterogeneity. Quartz-rich samples display major peaks at less than or equal to a 10 nm pore size, with a smaller effect on the PSD heterogeneity in most cases. In other brittle mineral-rich samples, there are no obvious major peaks, and a weak heterogeneity of the PSDs is displayed. (3) A greater TOC content, maturity, clay content and pore size can cause stronger heterogeneity of the PSD and higher fractal dimensions in the shale samples. This study helps to understand and compare the PSD and fractal characteristics from different samples and provides important theoretical guidance and a scientific basis for the exploration and development of shale gas resources.