The Poisson-Nernst-Planck (PNP) model plays an important role in simulating nanopore systems. In nanopore simulations, the large-size nanopore system and convection-domination Nernst-Planck (NP) equations will bring convergence difficulties and numerical instability problems. Therefore, we propose an improved finite element method (FEM) with an inverse averaging technique to solve the three-dimensional PNP model, named inverse averaging FEM (IAFEM). At first, the Slotboom variables are introduced aiming at transforming non-symmetric NP equations into self-adjoint second-order elliptic equations with exponentially behaved coefficients. Then, these exponential coefficients are approximated with their harmonic averages, which are calculated with an inverse averaging technique on every edge of each tetrahedral element in the grid. Our scheme shows good convergence when simulating single or porous nanopore systems. In addition, it is still stable when the NP equations are convection domination. Our method can also guarantee the conservation of computed currents well, which is the advantage that many stabilization schemes do not possess. Our numerical experiments on benchmark problems verify the accuracy and robustness of our scheme. The numerical results also show that the method performs better than the standard FEM when dealing with convection-domination problems. A successful simulation combined with realistic chemical experiments is also presented to illustrate that the IAFEM is still effective for three-dimensional interconnected nanopore systems.