In view of the problems of grouting pressure with multiple boreholes, according to Fourier series technique and correlative property of singular function, and based on natural boundary element method, the pressure boundary integral formula for Darcy flow equation of interior circular region with Dirichlet boundary value problem and the seepage pressure field of exterior vertical shaft with Neumann boundary value problem are deduced. Based on pressure boundary integral formula, the grouting pressure distribution of multi-hole grouting problems with arbitrary forms of hole-sitting are calculated. The numerical calculation shows that: the impermeability of the shaft wall produce large influence on the grouting pressure around the vertical shaft, accordingly, the shaft wall will suffer the bad influence. The pressure boundary integral formula of exterior vertical shaft will be more accurate compared with the Dirichlet problems by neglecting the impermeability effect.