Power to gas units and gas turbines have provided considerable opportunities for bidirectional interdependency between electric power and natural gas infrastructures. This paper proposes a series of multi-step strategy with surrogate Lagrange relaxation for operation co-optimization of an integrated power and natural gas system. At first, the value of coordination capacity is considered as a contract to avoid dysfunction in each system. Then, the uncertainties and risks analysis associated with wind speed, solar radiation, and load fluctuation are implemented by generating stochastic scenarios. Finally, before employing surrogate Lagrange relaxation, the non-linear and non-convex gas flow constraint is linearized by two-dimension piecewise linearization. In the proposed procedure, constraints for energy storages and renewable energy sources are included. Two case studies are employed to verify the effectiveness of the proposed method. The surrogate Lagrange relaxation approach with coordination branch & cut method enhances the accuracy of convergence and can effectively reduce the decision-making time.