I've given an answer below that is based on the WhiteboxTools geoprocessing library (for which I am the developer), but I would guess that you could use a similar workflow in just about any GIS. The only fairly unusual tool is the 'Depth-In-Sink' tool.
from whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"
in_dem = "DEM.tif"
# Invert the DEM
wbt.multiply(in_dem, -1, "iDEM.tif")
# Find the depressions, giving each one a unique ID
# Calcuate the depth-in-sink
# Assign each depression its maximum depth-in-sink value
wbt.extract_raster_statistics("depth_in_sink.tif", "depressions.tif", "max_dis.tif", stat="maximum")
# Threshold max_dis and create Booleans for deep deps vs everything else
wbt.greater_than("max_dis.tif", 50.0, "deep_depressions.tif", incl_equals=True)
wbt.less_than("max_dis.tif", 50.0, "not_deep_depressions.tif")
# Now fill the whole inverted DEM
wbt.fill_depressions("iDEM.tif", "iDEM_filled.tif", fix_flats=False)
# Multiply iDEM_filled by deep_depressions
wbt.multiply("iDEM_filled.tif", "deep_depressions.tif", "temp1.tif")
# Multiply iDEM by not_deep_depressions
wbt.multiply("iDEM.tif", "not_deep_depressions.tif", "temp2.tif")
# Add the two temp rasters together
wbt.add("temp1.tif", "temp2.tif", "temp3.tif")
# Invert the result back into the upward direction
wbt.multiply("temp3.tif", -1, "DEM_spikes_removed.tif")
If the features are true spikes (i.e. single grid cells with anomalously high values compared with neighbours), as is often the case with LiDAR DEMs, then you might consider using an image filtering approach instead of the depression-filling method. Filters that might be of use for this purpose (also available in WhiteboxTools) include the AdaptiveFilter, ConservativeSmoothingFilter, LeeFilter, MedianFilter, or the OlympicFilter. If your DEM is quite noisy, you might consider applying the FeaturePreservingDenoise tool. All of these tools may help to improve your DEM with a single operation compared to the above approach. If however your spikes are 'fat' (multiple cells in size) then the depression-filling approach may be best.
Lastly, you might consider instead filtering these points out of the LiDAR point cloud from which the DSM was created. The LidarRemoveOutliers tool will remove LiDAR points that exceed a specified height difference from the average elevation within a neighbourhood. This could be a very effective approach.