In this paper, the ground state Wigner function of a many-body system is explored theoretically and numerically. First, an eigenvalue problem for Wigner function is derived based on the energy operator of the system. The validity of finding the ground state through solving this eigenvalue problem is obtained by building a correspondence between its solution and the solution of stationary Schrödinger equation. Then, a numerical method is designed for solving proposed eigenvalue problem in one dimensional case, which can be briefly described by i) a simplified model is derived based on a quantum hydrodynamic model Cai et al. (2013) [7] to reduce the dimension of the problem, ii) an imaginary time propagation method is designed for solving the model, and numerical techniques such as solution reconstruction are proposed for the feasibility of the method. Results of several numerical experiments verify our method, in which the potential application of the method for large scale system is demonstrated by examples with density functional theory.