The integration of natural gas in electricity network requires a more reliable operating plan for increasing uncertainties in the whole system. In this paper, a threestage robust optimization model is proposed for resilient operation of energy system which integrates electricity and natural gas transmission networks with the objective of minimizing load curtailments caused by attacks. Nonconvex constrains are linearized in order to formulate the dual problem of optimal energy flow. Then, the proposed three-stage problem can be reformulated into a two-stage mixed integer linear program (MILP) and solved by Benders decomposition algorithm. Numerical studies on IEEE 30-bus power system with 7-node natural gas network and IEEE 118-bus power system with 14-node natural gas network validate the feasibility of the proposed model for improving resilience of integrated energy system. Energy storage facilities are also considered for the resiliency analysis.