dem - Interpolating LiDAR scattered ground data?

I wonder if someone can give me a hint on how to interpolate sparse data. I have a 1km by 1km tile located in a dense rain forest area.

Original LiDAR data are 37,395,021 points. I have played with them in several LiDAR processing programs (LAStools, MCC-LIDAR, Fusion, PDAL, to name a few) to get bare earth data.

Since the area covered has trees between 20-50m in height I have only 10% of the original data as ground points. For Instance, in one of my experiments I have 3,110,079 ground returns. The original data has a density of almost 27 points per square meter, while the ground data has 4 points per square meter. I have interpolated the data using SAGA Multilevel B-spline interpolation with a cell size of .5, 1 and 2m. I have also tried kriging, IDW, etc, and cannot obtain a smooth surface.

This is an example of the Multilevel B-spline interpolation showing ground points:

This is the same data using kriging:

Any ideas are welcome.