Your data is composed of a majoritarymajoriy of low-activity regions, with higher sparse-amplitiude sparse/concentrated activities. The bad news is that most of the methods to achieve your goal wouldare likely to somehow filter your data. The good news is that non-linear filtering could help you.
If you are not concerned with calculations and real-time yet, I would suggest first a centered median filter to extract a robust $0 g$ estimate, and remove it from the data. In other words: choose an appropriate left-right span of $K$. In each window around index $k$, computes the signal $\hat{s}(k) = s(k) - \textrm{median} [s_{k-K},\ldots,s_{k+K}]$.
A median is a non-linear filter minimizing the $ell_1$ norm (like the $ell_2$ norm minimization yields the standard mean). It is somewhat robust to outliers, and since yours are both above and below your thought $0 g$ estimate, a running median (computed around each sample) could extract this 0-level reference.
Looking at your plots, start with $K=50,100,200$ to see if this starts to provides you with what you are looking for. Then, more involded techniques could be considered.