Objective: With the growing complexity of modern power networks, coupled with the increasing frequency and intensity of natural hazards, there is an urgent need to develop advanced resilience assessment methods. Unlike conventional single-stage disaster events, multi-stage hazards such as earthquake-landslide-debris flow chains have successive and compounding impacts on electrical infrastructure. Such cascading events introduce severe uncertainty into the damage mechanisms and recovery processes of power systems, thereby threatening the continuity of electricity supply in disaster-affected regions. Traditional simulation-based approaches, especially those relying heavily on Monte Carlo techniques, often fail to capture the full dynamics of multi-stage shocks because of limited sampling or become computationally prohibitive when scaled to larger networks. To overcome these challenges, this study developed a probabilistic resilience assessment model that combines analytical probability calculations with optimization-based network analysis. The aim was to develop an accurate and efficient method for evaluating resilience in power systems under multi-stage disaster shocks. Methods: The proposed framework incorporates several method ological innovations. First, the entire power network is partitioned into smaller units through a community detection strategy, ensuring computational tractability while retaining the integrity of network interdependencies. Within these units, a mixed-integer programming model generates feasible power supply paths for each load node under normal operating constraints. This optimization-based representation identifies not only primary supply routes but also redundant paths that become critical in the event of failures. Second, the model incorporates probabilistic damage forecasting of transmission towers, which are among the most vulnerable components during seismic and secondary hazards. Instead of relying on scenario sampling, the model derives closed-form probability distributions of interruption durations during disasters and post-disaster recovery times. These distributions are obtained by analytically linking tower damage probabilities with repair processes, assuming realistic restoration practices. Finally, the resilience indicator is defined as the expected cumulative load-serving capability over the entire disaster cycle. By integrating the temporal evolution of service continuity, the framework captures both the degradation during hazard propagation and the recovery trajectory once repair efforts commence. This analytic approach eliminates the need for extensive sampling and significantly accelerates resilience estimation. Results: The proposed framework was validated using a modified IEEE 123-bus power network as the research object. This grid comprised one generation unit, 123 load centers, and 125 transmission lines. The application of the proposed method to this grid yielded several key research findings: First, the resilience of the power network exhibited an overall downward trend with increasing mainshock magnitude. It declined slowly under low-magnitude earthquakes with favorable disaster resistance. As the mainshock magnitude intensified and more chain-disasters occurred, the failure probability of transmission towers and lines increased, power supply redundancy was gradually weakened, and the network resilience dropped rapidly. Second, the established probabilistic resilience assessment model achieved markedly higher computational efficiency while guaranteeing satisfactory accuracy. By adopting analytical calculations, it realized accurate and efficient evaluation under the spatiotemporal evolution of earthquake disasters, effectively overcoming the excessive computational burden of Monte Carlo simulation. The resilience assessment results also revealed important system-level patterns: The resilience index of the studied grid was approximately 88.644 7% under a magnitude 7.0 earthquake, decreased to 76.854 9% under a magnitude 7.5, and further decreased to 58.413 4% under a magnitude 8.0. Although the resilience index declined significantly with increasing earthquake magnitude, the grid's supply path redundancy capability ensured that most electricity demand could still be met even under severe seismic events. These findings confirmed the crucial role of supply path redundancy and transmission tower reliability in maintaining stable electricity supply under cascading hazards. Conclusions: This study develops a novel probabilistic resilience assessment framework tailored for power networks exposed to multi-stage disaster shocks. By combining power unit partitioning, supply path optimization, and analytical probability modeling, the method addresses the limitations of simulation-heavy approaches and ensures the accuracy and efficiency of resilience estimation. The case study demonstrates the capability of the proposed method to quantify the interplay between disaster evolution, infrastructure vulnerabilities, and system recovery, offering highly relevant insights for planning resilient electricity infrastructure. This approach is particularly valuable for decision-makers and emergency planners, as it supports rapid assessment without sacrificing precision. In conclusion, the proposed method represents a significant advancement in resilience modeling for power systems. This study not only validates the feasibility of analytical probability-based approaches but also sets the stage for further research on integrating adaptive recovery strategies and resource-constrained repair models. In future work, the approach may be extended to consider real-time data integration, simultaneous restoration during hazard evolution, and multi-resource coordination, thereby enhancing the practical applicability of resilience assessment in real-world emergency contexts.