I don't think there's general formula or simplification for that. You might be able to find an approximate solution using perurbation techniquest if kernel $f(x)$ is very small. Or you can try to use regularization techniques such as Tikhonov regularization or total variation regularization to address the artifacts from deconvolution.

Thanks! Do you have any experience using perturbation to solve things like that? With perturbation $f = \delta + \lambda$, I'm still stuck with a $\int_{-\infty}^{\infty} g (\lambda * g)dx$ term...

I've tried myself and it seems that it's not sufficient to recover directly using perturbation. I think regularization might be best option for you

Deleted a comment based on misreading your question.

If you square the convolution it will be convolution in the Fourier domain. In that case the DC component will be a convolution and not just a multiplication of DC components. Maybe you can rewrite the squaring of the convolution as a convolution.

@mathreadler This transformation is not a simple multiplication of the Fourier transforms of the signals involved. Squaring a convolution is equivalent to the convolution of the autocorrelations of the signals. And it's still won't help in recovering $\int g^2 dx$ because after rewriting expression in terms of the Fourier transforms it will still involve both original and the conjugate Fourier transforms of the signals.

Is there any way to argue that $\int\, f(x-y)f(x-z)\, dx$ is proportional to $\delta(y-z)$?

@ruralreader it's true for some cases, e.g. $f(x)$ is Gaussian function

I can assume that $f(x)$ is Gaussian if it simplifies the problem enough. Usually $f(x)$ is either a 2D Gaussian or an Airy disk, but I'd be okay with whatever residuals I get assuming Airy is approximately Gaussian. (Sorry for the sudden introduction of 2D, I'd assumed the problem would be similar enough in 1 and 2 dimensions)

@LizaSazonova does the std $\sigma_f$ known?

@rumathe yup, approximately

What if we try to relate the integral of the convolution squared to the squared Fourier transforms of the functions being convolved then? Look's like Parseval's theorem.

But this will only work for gaussian $f(x)$ with known (approximate) std.

Using Fourier transform I can get to $\int (f*g)^2 dx = \big[ \int \hat{f} \hat{g} dk \big]^2$, is there a property of a Gaussian $f(x)$ we can use to simplify that?

What about @ruralreader's comment, is it true that $\int f(x-y)f(x-z)dx = \delta(y-z)$ if $f$ is Gaussian? In that case, we can re-write the convolution as $$ \int \int \int dt_1 dt_2 dx g(t_1) g(t_2) f(x-t_1) f(x-t_2) = \int \int dt_1 dt_2 g(t_1) g(t_2) \delta(t_1 - t_2)$$ which is easy to evaluate

well, yes @ruralreader's comment is correct. $$\int f(x-y)f(x-z)dx = \frac{1}{2\pi\sigma_f^2}\int \exp\bigg(-\frac{u^2+(y-z)^2/4}{2\sigma_f^2}\bigg)du = \frac{1}{\sqrt{2\pi}\sigma_f}\exp\bigg(-\frac{(y-z)^2}{4\sigma_f^2}\bigg) = \delta(y-z)$$

The problem here is that this property only apply to the result of convolution, not to the product of Gaussian functions.

Are you trying to restore true image of the object from e.g. telescope?

@LizaSazonova: Your line of reasoning and the other sundry comments lead to, for a tight low-pass filter (small $\sigma$), the result $\int g^2$ as an approximation.