Immersed boundary (IB) method has been extensively researched and used due to its advantage of simulating flow over complex and moving objects in rivers. However, for high Reynolds number flows, it is still a challenge to produce correct and smooth wall shear stress that is essential for sediment transport. This challenge comes from the non-uniform wall distance of IB cell center and the discontinuous wall functions in boundary layers. Specifically, when the background grid is fixed and the immersed boundary geometry is arbitrary, the distance from IB cell center to immersed boundary is non-uniform. Also, the distance to its image point, which is generally proportional to the distance of IB cell center, has the same non-uniformity and the image point may be located in the log-law layer or laminar layer. However, in many existing IB method, there are two discontinuous wall functions in these two boundary layers. As a result, the shear velocity may be estimated with image points and different wall functions, which gives a non-smooth wall shear distribution. Moreover, some of these distances are extremely small, resulting in numerical instability and even divergence. In this work, an improved treatment is proposed where the image point distance is set to be proportional to the size of IB cell, which produces a more uniform and stable distribution of the distance. Consequently, all image points are able to be in the log-law layer and follow the same wall function to provide smooth wall shear stress. This IB method has been implemented in the U2 RANS code. Cases in 1D and 2D are shown.