Perhaps you will find this a useful answer, at least for the symmetric case: Let $S^2_r(TM)\subset S^2(TM)$ denote the subbundle consisting of the cometrics of rank $r$. (This bundle is the disjoint union of $r{+}1$ smooth subbundles that are distinguished by the algebraic type of the symmetric bivector. However, this 'type division' won't play a role in this analysis.)
There is a natural first-order differential operator
$$
\delta_r:C^\infty\bigl(S^2_r(TM)\bigr)\longrightarrow
C^\infty\bigl(\Lambda^2(TM)\otimes \Lambda^r(TM)\otimes \Lambda^{r+1}(TM)\bigr)
$$
with the property that $\delta_r(g)=0$ if and only if the $r$-rank subbundle $g^\sharp(T^\ast M)\subset TM$ is integrable.
This differential operator is defined as follows: Locally, $g$ can be written in the form
$$
g = \sum_{i,j=1}^rg^{ij}\ X_iX_j\ ,
$$
where $X_1,\ldots, X_r$ are linearly independent vector fields, $g^{ij}=g^{ji}$ are functions, and $\Delta = \det (g^{ij})$ is nonvanishing. Now define (using the summation convention)
$$
\delta_r(g) = \Delta\ g^{ik}g^{jl}\ \left(X_i{\wedge}X_j\ \otimes\ X_1{\wedge}X_2{\wedge}\cdots{\wedge}X_r
\ \otimes\ [X_k,X_l]{\wedge}X_1{\wedge}X_2{\wedge}\cdots{\wedge}X_r\right).
$$
It is easy to verify that this is independent of the choice of basis $X_i$, and so it is globally defined. It is also obvious that it vanishes if and only if the bundle $g^\sharp(T^\ast M)$ (which is spanned locally by the $X_i$) is integrable.
Note that $\delta_r$ vanishes identically unless $2 \leq r < n$, as would be expected. Also, $\delta_r(\lambda g) = \lambda^{r+2}\ \delta_r(g)$ for all nonvanishing functions $\lambda$.
If I had considered oriented cometrics of rank $r$ (i.e., I had, in addition, fixed an orientation of $g^\sharp(T^\ast M)$), then I could then have got rid of the $\Lambda^r(TM)$ factor in the definition of $\delta_r$. However, I didn't want to impose the orientability condition, and the $\Lambda^r(TM)$ factor fixes that problem.
The almost Poisson case: Although you already have a criterion in the almost Poisson case, note that there is a similar construction to test the integrability of the bundle $\Pi^\sharp(T^\ast M)$ in the case of an almost Poisson structure $\Pi$ of constant half-rank $r$: Let $\Lambda^2_r(TM)\subset \Lambda^2(TM)$ be the subbundle of bivectors of half-rank $r$, i.e., its sections are those $\Pi$ such that $\Pi^r$ is nonvanishing as a section of $\Lambda^{2r}(TM)$ while $\Pi^{r+1}\equiv0$.
There is a natural first-order differential operator
$$
\partial_r:C^\infty\bigl(\Lambda^2_r(TM)\bigr)\longrightarrow
C^\infty\bigl(\Lambda^2(TM)\otimes\Lambda^{2r+1}(TM)\bigr)
$$
defined as follows: Locally, write $\Pi$ in the form
$$
\Pi = \sum_{i,j=1}^{2r} {\tfrac12}\ a^{ij}\ X_i\wedge X_j\ ,
$$
where $X_1,\ldots,X_{2r}$ are linearly independent vector fields, $a^{ij}=-a^{ji}$ are functions, and $\mathrm{Pf}(a^{ij})$ is nonvanishing (since $\Pi^r = r!\ \mathrm{Pf}(a^{ij})\ X_1\wedge X_2\wedge\cdots\wedge X_{2r}\not=0$). Now define (again using the summation convention)
$$
\partial_r(\Pi) =
{\tfrac14}\ a^{ik}a^{jl}\
X_k{\wedge}X_l\otimes\ [X_i, X_j]{\wedge}\Pi^r\ .
$$
It is easy to verify that this does not depend on the choice of basis $X_i$,
and so it is well-defined.
It is evident that that $\partial_r(\Pi)$ vanishes if and only if the bundle $\Pi^\sharp(T^\ast M)$ (which is spanned locally by the $X_i$) is integrable. Also, $\partial_r(\lambda\ \Pi)= \lambda^{r+2}\partial_r(\Pi)$ for all nonvanishing functions $\lambda$.
This is a somewhat more explicit test that the one you proposed, I think.
The general case: It turns out that one does not need any assumption of symmetry or skew-symmetry to get the integrability tensor, just an assumption of constant rank. Let $\otimes^2_r(TM)\subset TM\otimes TM$ denote the subbundle consisting of the tensors of rank $r$. A section $\phi\in C^\infty\bigl(\otimes^2_r(TM)\bigr)$ can be written locally in the form
$$
\phi = \sum_{i,j=1}^r f^{ij}\ X_i\otimes Y_j
$$
where $X_1,\ldots,X_r$ are linearly independent vector fields, $Y_1,\ldots,Y_r$ are linearly independent vector fields, and $f=(f^{ij})$ is an invertible matrix, i.e., $\det(f)\not=0$. The vector fields $X_i$ are local sections of a globally defined bundle $\lambda_\phi\subset TM$ of rank $r$, and the vector fields $Y_i$ are local sections of a globally defined bundle $\rho_\phi\subset TM$ of rank $r$.
There is a natural first-order differential operator
$$
D_r: C^\infty\bigl(\otimes^2_r(TM)\bigr)\longrightarrow
C^\infty\bigl(\Lambda^2(TM)\otimes \Lambda^r(TM)\otimes \Lambda^{r+1}(TM)\bigr)
$$
which, relative to the local expression given above, takes the form
$$
D_r(\phi) =
\det(f)\ f^{ik}f^{jl}\ \left(Y_k{\wedge}Y_l\ \otimes
\ Y_1{\wedge}Y_2{\wedge}\cdots{\wedge}Y_r
\ \otimes\ [X_i,X_j]{\wedge}X_1{\wedge}X_2{\wedge}\cdots{\wedge}X_r\right).
$$
This operator has the property that $D_r(\phi)=0$ if and only if $\lambda_\phi$ is integrable.
Note that there is a well-defined involution $\iota$ on $C^\infty\bigl(\otimes^2_r(TM)\bigr)$ such that
$$
\iota(\phi) = \sum_{i,j=1}^r f^{ij}\ Y_j\otimes X_i\ ,
$$
and the integrability of $\rho_\phi$ is equivalent to the equation $D_r\bigl(\iota(\phi)\bigr)=0$.
