This paper proposes a new structural topology optimization method using a dual-level point-wise density approximant and the meshless Galerkin weak-forms, totally based on a set of arbitrarily scattered field nodes to discretize the design domain. The moving least squares (MLS) method is used to construct shape functions with compactly supported weight functions, to achieve meshless approximations of system state equations. The MLS shape function with the zero-order consistency will degenerate to the well-known 'Shepard function', while the MLS shape function with the first-order consistency refers to the widely studied 'MLS shape function'. The Shepard function is then applied to construct a physically meaningful dual-level density approximant, because of its non-negative and range-restricted properties. First, in terms of the original set of nodal density variables, this study develops a nonlocal nodal density approximant with enhanced smoothness by incorporating the Shepard function into the problem formulation. The density at any node can be evaluated according to the density variables located inside the influence domain of the current node. Second, in the numerical implementation, we present a point-wise density interpolant via the Shepard function method. The density of any computational point is determined by the surrounding nodal densities within the influence domain of the concerned point. According to a set of generic design variables scattered at field nodes, an alternative solid isotropic material with penalization model is thus established through the proposed dual-level density approximant. The Lagrangian multiplier method is included to enforce the essential boundary conditions because of the lack of the Kronecker delta function property of MLS meshless shape functions. Two benchmark numerical examples are employed to demonstrate the effectiveness of the proposed method, in particular its applicability in eliminating numerical instabilities. when the prescribed objective function is optimized under design constraints. So far several typical methods have been developed for topology optimization, for example, the homogenization method [2], solid isotropic material with penalization (SIMP) approach [3,4], and level set-based method [e.g., 5-8].However, topology optimization of continuum structures essentially belongs to a family of integer programming problems with a large number of discrete design variables (0,1), to which many more efficient gradient-based optimization algorithms cannot be directly applied [1]. To this end, the homogenization and SIMP methods are the widely used approaches to relax the original problem, to allow discrete design variables to take intermediate density values ranging from 0 and 1. In particular, SIMP, as an extension of the homogenization method, has received much popularity in the community of structural optimization because of its conceptual simplicity and implementation easiness. With SIMP, the original (0,1) discrete optimization is changed to the optimization with ...