With the growing maturity of multi-robot system technology, its applications have expanded across various domains. This paper addresses the critical issue of task allocation in nuclear accident rescue scenario, which plays a pivotal role in the success of such operations. The problem is formulated as a multi-objective optimization problem, taking into account three key indicators: execution time, radiation accumulation, and waiting cost. To effectively tackle this problem, an two-stage evolutionary algorithm is proposed. Firstly, a solution encoding method and a crossover mutation method is devised tailored to the problem’s characteristics. Secondly, a two-stage search strategy is designed. In the first stage, a fixed population size and shift-based density estimation method are used to quickly converge the solution set to the Pareto front. The latter stage uses an infinite size population to find as many Pareto solutions as possible. Finally, a local search strategy is introduced to improve the quality of solution set. In the experimental section, our proposed method is compared with five state-of-the-art algorithms on nine instances of varying scales. Across five evaluation metrics, the proposed algorithm demonstrates competitive performance on all instances. These results underscore the efficacy and competitiveness of our approach in tackling the task allocation problem in multi-robot systems within nuclear accident rescue.