One other suitable open source tool is GMT. The "surface" algorithm makes pretty nice results and it supports so called soft breaklines https://docs.generic-mapping-tools.org/latest/surface.html.
Use x, y, z data in the breakline file as a ‘soft breakline’. A ‘soft
breakline’ is a line whose vertices will be used to constrain the
nearest grid nodes without any further interpolation. A coastline or a
lake shore are good examples of ‘soft breaklines’. Multi-segments
files are accepted. If your lines do not have z-values or you wish to
override those with a constant z-value, then append +zlevel to the
filename. If no value is given then we default to 0.
As x, y, z data for example POLYGONZ type shapefiles which have 0 as z value can be used. The hardest thing with GMT is to learn the command line syntax.
The GRASS function v.surf.rst is quite similar to the GMT surface https://grass.osgeo.org/grass82/manuals/v.surf.rst.html. This function needs the shoreline as a raster mask:
User can either use r.mask to set a mask or specify a raster map in
mask option, which will be used as a mask. The approximation is
skipped for cells which have zero or NULL value in mask. NULL values
will be assigned to these cells in all output raster maps. Data points
are checked for identical points and points that are closer to each
other than the given dmin are removed. If sparsely digitized contours
or isolines are used as input, additional points are computed between
each 2 points on a line if the distance between them is greater than
specified dmax. Parameter zmult allows user to rescale the values used
for approximation (useful e.g. for transformation of elevations given
in feet to meters, so that the proper values of slopes and curvatures
can be computed).
You can also try to extract the vertices of the shoreline as points with 0 depth and append them into your sonar data. Then you can use your current interpolation. The drawback compared to breakline or mask is that the interpolation algorithm does not know that the shoreline points are somehow special and therefore the algorithms tend to overshoots the interpolation over the dry land, and zero depths are not necessarily kept as zeroes.