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](https://mathoverflow.net/questions/290220/limit-associated-with-complementary-sequences) or [here](https://mathoverflow.net/questions/290808/a-possibly-surprising-appearance-of-sqrt2?noredirect=1&lq=1), 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 less clear to me, but it may possibly be studied in the same lines.