In this paper, we investigate the coupling of natural boundary element and finite element methods of exterior initial boundary value problems for hyperbolic equations. The governing equation is first discretized in time, leading to a time-step scheme, where an exterior elliptic problem has to be solved in each time step. Second, a circular artificial boundary $\Gamma_R$ consisting of a circle of radius $R$ is introduced, the original problem in an unbounded domain is transformed into the nonlocal boundary value problem in a bounded subdomain. And the natural integral equation and the Poisson integral formula are obtained in the infinite domain $\Omega_2$ outside circle of radius $R$. The coupled variational formulation is given. Only the function itself, not its normal derivative at artificial boundary $\Gamma_R$, appears in the variational equation, so that the unknown numbers are reduced and the boundary element stiffness matrix has a few different elements. Such a coupled method is superior to the one based on direct boundary element method. This paper discusses finite element discretization for variational problem and its corresponding numerical technique, and the convergence for the numerical solutions. Finally, the numerical example is presented to illustrate feasibility and efficiency of this method.