This is less more than a long comment. It seems that a somehow simpler case is the when $a_{n+1}-a_n>1$ (at least eventually) which means $b_{n+1}-b_n\in\{1,2\}$ (eventually), and $a_{n+1}-a_n$ takes values in a finite set $\mathcal{A}$ of say $r$ integers greater or equal to $2$. Then the sequence $a_{n+1}-a_n$ is generated by a substitution map (like in the one described here or here, so that the entire sequence $a_{n+1}-a_n$, seen as an infinite string, has the form $p\cdot\tau(p)\cdot\tau(p)^2\cdot\tau(p)^3\dots$, where the dots stand for concatenation of words, and $p$ is a prefix. For such a substitution map $\tau$ one can write an associated $r\times r$ transition matrix $A:=(a_{ij})$ for the number of occurrences of each symbol in a transformed word: let $a_{ij}$ denote the number of occurrences of the letter $s_i$ in the word $\tau(s_j)$. Then, if in the vector $X$, $X_i$ is the number of occurrences of $s_i$ in a given word $w$, the coordinates of $AX$ give the number of occurrences of each $s_i$ in the word $\tau(w)$, so that, for $w=p$, the vector $A^k X$ gives the distribution of letters in $\tau(p^k)$; the length of $\tau(p^k)$ is the sum of coordinates of $A^{k-1} X$. This way the various asymptotics of $a_n$ can be easily related to the spectrum of $A$; in fact one use the standard techniques of finite Markov chains introducing a suitable Markov chain deduced from the map $\tau$.
The case when $a_{n+1}-a_n$ can take the value $1$ frequently, so that $b_{n+1}-b_n$ assumes more values, seems more clear to me, but it may possibly be studied in the same lines.