According to Peskin & Schroeder (page 325), the Feynman rule for the counterterm
------(x)-----
for
$$ \frac12 \delta_Z(\partial_\mu\phi_r)^2-\frac12\delta_m \phi_r^2$$
being $\phi_r$ the renormalized field, is given by
$$i(p^2\delta_Z-\delta_m)$$
which resembles rather the (multiplicative) inverse of the propagator for the original Lagrangian (whith physical quantities). Why?
Counter Term Feynman Rules
A counter term in a Lagrangian is interpreted as an interaction, and the Feynman rules may be derived by the Feynman path integral in the standard way. This accounts for the factor of $i$ and the functional derivatives will take care of the factor of $\frac12$.
Just like in $\frac{\lambda}{4!}\phi^4$ theory we have a $-i\lambda$ Feynman rule, with $\frac{\delta_m}{2}\phi^2$ we have a $-i\delta_m$ Feynman rule. In addition, derivatives of fields give rise to momenta: $(\partial \phi)^2 \mapsto p^2$. So the overall counter term is,
$$i(p^2\delta_Z-\delta_m).$$
The key is, if in doubt, always check things explicitly by treating the terms as interactions in the Feynman path integral.
Alternative Method
I highly recommend Collins' book on renormalisation for a method of renormalisation which dispenses with a Lagrangian and is ideal for systematically computing diagrams. This involves a combination of the BPHZ method and dimensional regularisation.
Essentially, the Zimmermann forest formula generates all the counter term graphs. Then instead of inserting Feynman rules, you insert a particular 'subtraction' operator to remove the divergences.
Of course, this can be related back to a Lagrangian by considering what the Feynman rules would read, and equating the two results to find the counter terms to whatever order you need.
Thee diverging terms for the propagator come from the renormalization of the self-energy $\Sigma$, defined by $G^{-1}=G_0^{-1}-\Sigma$, where $G_0$ is the propagator defined by the Lagrangian (i.e. bare propagator + counterterms) :
$G^{-1}_0=(1+\delta Z)p^2+(m^2_0+\delta m^2)$.
One chooses the counterterms to cancel the divergences coming from $\Sigma$ order by order. If the theory is perturbatively renormalizable, only these two counterterms are sufficient at every order in perturbation theory.
I also found it hard to find a rigorous explanation. By dimensional analysis and explicit checking you indeed find that this is the correct Feynman rule for the counterterm.
The closest thing I have to a derivation is that we can include the counterterm in the kinetic term and expand this full propagator for small $\delta Z$
$ G = \frac{i}{(1+\delta Z)p^2+m^2} \approx \frac{1}{p^2+m^2}-\frac{\delta Z}{(p^2+m^2)^2} + ...$.
I think that not including the $(p^2+m^2)^2$ term in the Feynman rule is correct since this counterterm gets treated as an (amputated) vertex. A four point function similarly does not include 4 propagator factors in its feynman rule.
I do not know why I get a different sign than Peskin & Schroeder.
Also, the expansion is questionable because $\delta Z$ is in fact infinite as opposed to small.
