The regular hybrid boundary node method is improved and a new numerical method in geotechnical engineering is proposed in this paper. The source points of the fundamental solution are located on the boundary directly. There is no free parameter in the method which combines the modified variational function with the moving least squares approximation. It not only has the well-known dimension-reduction advantage of the boundary element method

but also is a true meshless method. Only randomly distributed nodal points on the bounding surface of the domain are required and the post-processing is very simple. At the same time

it is suitable for treatment of the semi-infinite space problems in geotechnical engineering. The numerical examples show that this method has high precision and is more suitable for computation in geotechnical engineering.