We develop a parallel two-level domain decomposition method for the 3D unsteady source identification problem governed by a parabolic partial differential equation (PDE). The domain of the PDE is firstly decomposed into several overlapping subdomains and the original inverse source identification problem is then transformed into smaller independent subproblems defined on these subdomains. Each subproblem is formulated as a PDE-constrained optimization problem with appropriate conditions prescribed on the inner boundaries and discretized by finite element method. The resulting coupled algebraic systems are solved simultaneously by restarted GMRES method with a space-time restricted additive Schwarz preconditioner. When forming the preconditioner, a second level of domain decomposition is introduced for each subdomain. The solutions of these subproblems are combined together to form an approximated global solution to the original inverse problem by discarding the overlapping parts of the solution. Since all the subproblems are solved independently, the two-level domain decomposition method provides higher degree of parallelism and saves much computing time. Numerical experiments conducted on a supercomputer with thousands of processor cores validate the efficiency and robustness of the proposed approach.