Natural fractures in middle–deep shale reservoirs experience continuous stress evolution under the combined influence of hydraulic stimulation and long-term production. Their stability evolution critically affects wellbore integrity, stimulation effectiveness, and the safety of infill-well deployment. Focusing exclusively on natural fractures, this study develops a stability evaluation framework by coupling a 3D geomechanical model with a discrete fracture network (DFN). The approach maps the strike and dip of natural fractures onto the 3D grid and constructs a spatially continuous fracture-orientation volume through geostatistical interpolation. This enables the unified coupling of natural fracture geometry with the regional 3D stress field and rock mechanical attributes, providing a continuous 3D quantification of natural-fracture stability. Results show that natural-fracture slip risk is governed by the combined effects of fracture orientation and tectonic stress regime, with distinct high-risk orientations under normal-faulting, strike-slip, and reverse-faulting conditions. Fluid injection may trigger natural-fracture instability through reduced effective normal stress and lowered frictional strength, whereas long-term production enhances effective stress and generally improves natural-fracture stability. In the YS108 block, the stress regime evolves toward a typical normal-faulting state after nearly ten years of production, leading to significantly reduced slip risk of natural fractures. The proposed 3D evaluation framework provides a practical basis for post-production stability assessment, well-trajectory optimization, and stimulation-risk management in middle–deep shale reservoirs.