The calculation of inter-granule contact force in three-dimensional (3D) granular systems is a key and challenging aspect of granular mechanics research. Two elastic rubber balls are used as research objects for in-situ flat pressing micro-CT experiments. Based on the Hertzian contact theory and Tatara large deformation contact theory, the contact model of elastic balls is verified, and the theoretical formula of the contact force of elastic balls based on the experiment is obtained. Taking the 3D granular systems as research object, in-situ probe loading experiment of micro-CT is carried out to obtain the 2D image sequence of the granules, after a series of digital transformations, the digital body images emerge, the contact force networks of the 3D granular systems under different loading conditions are obtained by constructing pore network models. The contact force distribution and evolution law of the granular systems are analyzed. The relation among the number of strong contacts, the distribution evolution, and the stability of the granular system is explored. The results show that the two elastic ball contact model conforms to the Hertzian contact theory and Tatara large deformation contact theory, and the contact force fitting formula based on experiment can characterize the contact force between two granules reasonably and effectively. The contact force of granules under probe loading is distributed in a net-like pattern starting from the contact point of the indenter and gradually transmitted to the lower and the surrounding area. The trend of average contact force is consistent with the trend of the contact times, showing a significant phase transition. With the increase of contact times, the frequency of particle compression increases, resulting in a greater contact force between granules, ultimately stabilizing at about 10.5 N. The number of strong contacts accounts for 45% to 50% of the total number of contacts, distributed throughout the whole granular system and supporting the network structure of the granular system. The larger values are concentrated below the indenter and exhibit a branching distribution. In the loading process, an equilibrium point is established at z = 14 mm, where the number of strong contacts reaches the peak. The network structure of strong contact force is spread throughout the entire 3D granular system, establishing the main skeleton that can withstand external loads. As the loading continues, the total value of strong contact forces increases, and their distribution within the granular system becomes more uniform.