Okay, I understand Nate's comment now and am posting this elucidation of it as a community wiki answer.

Recall that an ordered set partition of $X$ into $j$ blocks is an ordered tuple $(T_1,T_2,\ldots,T_j)$ of non-empty subsets $\varnothing \neq T_i \subseteq X$ that are pairwise disjoint and whose union is all of $X$. Since $G$ acts on $X$ it acts naturally on the ordered set partitions of $X$. Moreover, to any ordered set partition $(T_1,T_2,\ldots,T_j)$ and choice of subset $\{y_1 < y_2 < \cdots < y_j\} \subseteq Y$ of $j$ colors, we can associate the coloring $f\colon X \to Y$ determined by $f(x) = y_i$ iff $x \in T_i$. And two such colorings are equivalent under the action of $G$ exactly when the ordered set partitions are equivalent.

This shows $$ P(k) = \sum_{j\geq 1} (\textrm{$\#$ of orbits of $G$ acting on ordered set partitions of $X$ into $j$ blocks}) \cdot \binom{k}{j},$$ i.e., that $$ a_j = \textrm{$\#$ of orbits of $G$ acting on ordered set partitions of $X$ into $j$ blocks}.$$

This is a perfectly fine formula that shows the $a_j$ are nonnegative integers, but it's a bit different from the Burnside's Lemma formula for $P(k)$ in that the Burnside's Lemma formula is "local" in the sense of looking at each $g\in G$ individually. I guess it would be hard to give a similar "local" formula for the $a_j$.