A radial basis function ghost cell immersed boundary method of simulating flows around arbitrary complex or multiple immersed boundaries is proposed in this paper. In this method, incompressible Navier-Stokes equations are discretized on fixed Cartesian staggered gridby the finite difference method. A fractional step method is used for time integration, together with third order Runge-Kutta scheme. A high-order TVD MUSCL (total variation diminishing monotonic upstream-centered scheme for conservation law) scheme is used to discretize convective terms. Two salient features are emphasized in the present study. First, boundary conditions at the immersed interface are enforced by a continuous ghost cell method to consider the influence of immersed boundary on the flow field. The immersed bodies are treated as virtual boundaries immersed in the flow. And Navier-Stokes equations are solved in the entire computation domain, including solid domain. Therefore, programming complexity is greatly reduced and the treatment of immersed boundaries is simplified. Second, a polynomial and radial basis function is introduced to implicitly represent and reconstruct arbitrary complex immersed boundaries. Iso-surface distance functions about interface geometries are fitted with some sampling points of body surfaces. It is flexible and robust. Moreover, the information about interface positions on the background grid can be easily identified by the signed distance functions. Based on our in-house developed immersed boundary method solver, typical test cases are simulated to validate the proposed method. The flows around a cylinder at Reynolds numbers of 40, 100 and 200 are first simulated and a grid resolution study is carried out. Good agreement is achieved by comparing with previous numerical results, which shows that this method is accurate and reliable. In the second case of flow around airfoil, the good agreement with previous study shows that the present method has the ability to simulate complex immersed boundary flow. In the last case of flow around array of thirteen cylinders, the ability of present method for multiple immersed boundaries is well proved. And hydrodynamic interaction among multiple bodies is briefly analyzed.