In microwave imaging, a great chunk of literature is devoted to regularization and its effect on the solution process and inversion results. One of the common methods for microwave imaging is the Contrast Source Inversion (CSI) method, which is essentially a gradient-based optimization.
In CSI, one would formulate the inverse problem as the optimization of the functional $\mathcal F(\chi)$:
$$
\mathcal F_n(\chi) = F_{1,n}(\chi)+\ldots+ F_{m,n}(\chi)
\tag{1}
$$
In (1), $F_{i,n}$ is the $i$th cost function (in particular for CSI, $m=2$) for the $n$th iteration. Now, there are several ways to perform a regularization on (1).
- additive regularization with a penalty function $F^R(\chi)$ and weighting parameter $\gamma$:
$$
\mathcal F_n(\chi) = F_{1,n}(\chi)+\ldots+ F_{m,n}(\chi)+\gamma^2F^R_n(\chi)
\tag{2}
$$
- multiplicative regularization:
$$
\mathcal F_n(\chi) = \big(F_{1,n}(\chi)+\ldots+ F_{m,n}(\chi)\big)F^R_n(\chi)
\tag{3}
$$
The problem with additive regularization is that looking for good $\gamma$ usually requires a lot of numerical experimentation and, possibly, prior information about the desired reconstruction of $\chi$.
Multiplicative regularization stems from the idea that $\gamma$ should depend on the $\mathcal F_{n-1}(\chi)$, the functional value at the previous iteration.
I suggest looking into
In this free-access paper, the authors present a multiplicative regularizer that fits into the gradient-optimization framework.
Now, regarding the usual suspects from the regularizers' family:
$$
F^R_n(\chi) = \frac{\int\limits_S ds(\vec{\rho}) \sqrt{\left|\nabla \chi(\vec{\rho})\right|^2+\delta_n^2}}{\int\limits_S ds(\vec{\rho}) \sqrt{\left|\nabla \chi_n(\vec{\rho})\right|^2+\delta_n^2}}
\tag{4}
$$
- $L^2$-norm of the total variation - prefers smooth profile (that's probably what you are looking for), but you will lose the distinct object edges. Although, one can employ some edge-preserving algorithm, - but that's a topic for another discussion.
$$
F^R_n(\chi) = \frac{\int\limits_S ds(\vec{\rho}) \left(\left|\nabla \chi(\vec{\rho})\right|^2+\delta_n^2\right)}{\int\limits_S ds(\vec{\rho}) \left(\left|\nabla \chi_n(\vec{\rho})\right|^2+\delta_n^2\right)}
\tag{5}
$$
In (4) and (5), $\rho$ is the 2-D position vector (wrt to which we perform spatial differentiation and integration), $S$ is your 2-D domain, and $\delta_n$ can be defined, as follows:
$$
\delta_n^2 = h^2\sum\limits_{j}F_{j,n}(\chi)
$$
By choosing $j$ one can control, which cost functions from (1) participate in regularization. If you have just one functional ($m=1$), you have a simpler formula. $h$ is the reciprocal of your mesh element size.
Now, choosing non-zero $\delta_n$ helps to restore differentiability of regularizers $F_n^R$ in (4) and (5).
In this answer, I simplified notation from the cited paper to only one inverse function $\chi$ (CSI uses two different ones) and 2-D. The paper (and especially references in it) should help in finding regularizers that work better for your particular inverse problem or are just simpler versions of them.