Objectives  For spherical harmonic synthesis of magnetic field at lunar irregular surface, traditional calculating algorithm has a low computing efficiency and the possible instability. Therefore, the radial Taylor's series expansion method is proposed.

  Methods  Firstly, the formulae of the traditional and newly presented methods are given. Then, we adopt a spherical harmonic model of the lunar magnetic field with maximum degree and order of 450, the numerical tests are performed in the lunar highland region with a drastic topographic relief. Finally, the global magnetic fields on lunar relief surface and on the reference sphere are calculated, respectively.

  Results  The results of numerical tests and practical application show that the calculating accuracy of our method mainly depends on the maximum order of the Taylor's series expansion, and if the observing surface changes more dramatically and the spherical harmonic model of the magnetic field has higher maximum degree and order, a higher maximum order of the Taylor's series expansion is required. Moreover, the average radius of the algorithm is suggested to be the area-weighted mean of the radius values of all calculating points. The computing time varies linearly with the maximum order of the Taylor's series expansion. If more calculating points and higher maximum degree and order of spherical harmonic model of the magnetic field, the computational efficiency will be improved more significantly by our proposed method. The differences between magnetic fields on lunar relief surface and on the reference sphere are very large, which suggests that the source depth of the internal magnetic field of the Moon is relatively shallow. Especially, the differences of the intensity of the magnetic field have a negative correlation with the topography variations, but there is no correlation between the radial component of the magnetic field and the topography variations. These indicate that the sources' magnetized directions of the internal magnetic field of the Moon are not radial.

  Conclusions  In brief, this study significantly certifies that our proposed method has a high computing efficiency and a very weakly instability.