Although this is 7 years later, I noticed that there is some gap in the literature in how to derive Feynman rules for such type of Lagrangians. Let me type up how to derive the vertex for the Lagrangian you mentioned in order to help people in the future who work on QFT since I think the other answers on forums are a bit too unclear for beginners and not so systematic in general. Calculating Feynman rules boils down to computing propagators, from that one can construct generating function and include interactions. The Feynman rules then follow by taking a bunch of functional derivatives.
So let us start with the Lagrangian:
\begin{equation}
\mathcal{L}_{\mathrm{int}} = \frac{gh}{3!} \phi^3 \partial^2 \phi.
\end{equation}
The functional approach of QFT can be used to derive the vertex Feynman rule. Let me try to give a self-contained answer/strategy to such questions. Suppose that the scalar field $\phi$ obeys the standard Klein-Gordon free field Lagrangian $\mathcal{L}_0 = \frac{1}{2} (\partial_\mu \phi)^2 - \frac{1}{2}m^2 \phi^2$. Then the propagator can be found by writing the free part of the action as:
\begin{equation}
S_0 = \int d^4 x \Big[\frac{1}{2} (\partial_\mu \phi)^2 - \frac{1}{2} m^2 \phi^2\Big] = \int d^4 x \Big[\frac{1}{2} \phi (-\partial^\mu \partial_\mu - m^2)\phi\Big],
\end{equation}
where I used integration by parts for the first term. This suggests that the propagator $D(x-y)$ can be found by:
\begin{equation}
-(\partial^\mu \partial_\mu + m^2)D(x-y) = i \delta(x-y).
\end{equation}
The $i$ in here is just a matter of convention. In Fourier space we see quite easily that:
\begin{equation}
D(k) = \frac{i}{k^2 - m^2}.
\end{equation}
Let $J(x)$ be some function of spacetime which we call the current associated to $\phi$. The generating function can be written as a path integral:
\begin{equation}
Z[J] = \int \mathcal{D}\phi e^{i \int d^4 x (\mathcal{L} + J(x)\phi(x))},
\end{equation}
where $\mathcal{L} = \mathcal{L}_0 + \mathcal{L}_{\mathrm{int}}$ and $\mathcal{D}\phi$ denotes the integration over all field configurations (for a gauge field for instance you would need to be careful about fixing a gauge in this expression but for scalar field there are no issues). In this language it is simple to illustrate that the free part of the generating function can be written in the form:
\begin{equation}
Z_0[J] = Z_0[0] e^{-\frac{1}{2} \int d^4 x d^4 y J(x)D(x-y)J(y)}.
\end{equation}
This can be done by starting from the definition of $Z[J]$, only consider the free Lagrangian and define the shifted field $\phi^\prime(x) = \phi(x) - i \int d^4 y J(y)D(x-y)$. The subscript $0$ indicates that we consider the free part only.
In order to find the expression for the 4-point vertex you wrote, one must compute the 4-point correlation function:
\begin{equation}
\langle \Omega|T \{\phi(x_1) \phi(x_2) \phi(x_3) \phi(x_4)\}|\Omega \rangle = \frac{\int d^4 x \phi(x_1)\phi(x_2) \phi(x_3) \phi(x_4) e^{i \int d^4 x \mathcal{L}}}{\int d^4 x e^{i \int d^4 x \mathcal{L}}}.
\end{equation}
It is simple to check that in the language of functional derivatives and the generating function $Z[J]$ one could write this as:
\begin{align}
&\langle \Omega|T \{\phi(x_1) \phi(x_2) \phi(x_3) \phi(x_4)\}|\Omega \rangle = \frac{1}{Z[0]}\Big(-i \frac{\delta}{\delta J(x_1)}\Big)\Big(-i \frac{\delta}{\delta J(x_2)}\Big) \\
&\Big(-i \frac{\delta}{\delta J(x_3)}\Big)\Big(-i \frac{\delta}{\delta J(x_4)}\Big)Z[J] \Big|_{J=0}.
\end{align}
Okay so now we have everything ready for the computation. The total generating function can be written as:
\begin{equation}
Z[J] = \int \mathcal{D}\phi e^{i \int d^4 x \mathcal{L}_{\mathrm{int}}} Z_0[J] = \int \mathcal{D}\phi e^{i \int d^4 x \frac{gh}{3!} (-i\frac{\delta}{\delta J(x)})^3 \partial^2 (-i \frac{\delta}{\delta J(x)})}Z_0[J].
\end{equation}
If $g,h$ are small then we can look at a perturbative expansion of this exponential. Focusing on the interacting part at leading-order we find that:
\begin{equation}
Z[J] \approx i \int d^4 x \frac{gh}{3!}(-i\frac{\delta}{\delta J(x)})^3 \partial^2 (-i \frac{\delta}{\delta J(x)})Z_0[J].
\end{equation}
Therefore the correlation function can be written as:
\begin{equation}
\langle \Omega|T \{\phi(x_1) \phi(x_2) \phi(x_3) \phi(x_4)\}|\Omega \rangle = \Big(-i \frac{\delta}{\delta J(x_1)}\Big)\Big(-i \frac{\delta}{\delta J(x_2)}\Big) \\
\Big(-i \frac{\delta}{\delta J(x_3)}\Big)\Big(-i \frac{\delta}{\delta J(x_4)}\Big)i \int d^4 x \frac{gh}{3!}(-i\frac{\delta}{\delta J(x)})^3 \partial^2 (-i \frac{\delta}{\delta J(x)}) e^{-\frac{1}{2} \int d^4 x d^4 y J(x)D(x-y) J(y)}|_{J=0}.
\end{equation}
Let us exploit notation in which we leave the integral over $x$ implicit and e.g. $D_{xy} = D(x-y)$, $J_y = J(y)$, $D_{xy} J_y \equiv \int d^4 y D(x-y)J(y)$ to make expressions more concise. We will also leave the evaluation at $J=0$ implicit. Bubble diagrams, i.e. ones which include $D_{xx}$ will be ignored in the expansion below. Also we will ignore terms which will give zero for $J=0$.
So let us start computing this correlation function using this compact notation:
\begin{align}
&-i \langle \Omega|T \{\phi(x_1) \phi(x_2) \phi(x_3) \phi(x_4)\}|\Omega \rangle = \frac{gh}{3!} \frac{\delta}{\delta J_1} \frac{\delta}{\delta J_2} \frac{\delta}{\delta J_3}\frac{\delta}{\delta J_4} \Big(\frac{\delta}{\delta J_x}\Big)^3 \partial^2 \Big(\frac{\delta}{\delta J_x}\Big) e^{-\frac{1}{2}J_z D_{zy}J_y} \\
&= \frac{gh}{3!} \frac{\delta}{\delta J_1} \frac{\delta}{\delta J_2} \frac{\delta}{\delta J_3}\frac{\delta}{\delta J_4} \Big(\frac{\delta}{\delta J_x}\Big)^3 e^{-\frac{1}{2}J_z D_{zy}J_y} \partial^2 (-D_{xy}J_y) \\
&= \frac{gh}{3!} \frac{\delta}{\delta J_1} \frac{\delta}{\delta J_2} \frac{\delta}{\delta J_3}\frac{\delta}{\delta J_4} \Big(\frac{\delta}{\delta J_x}\Big)^2 J_z D_{xz} e^{-\frac{1}{2} J_z D_{zy}J_y} \partial^2 (D_{xy}J_y) \\
&= \frac{gh}{3!} \frac{\delta}{\delta J_1} \frac{\delta}{\delta J_2} \frac{\delta}{\delta J_3}\frac{\delta}{\delta J_4} \frac{\delta}{\delta J_x} \Big(-D_{xz} J_z D_{xv} J_v e^{-\frac{1}{2}J_z D_{zy}J_y} \partial^2 (D_{xy} J_y)\Big) \\
&= \frac{gh}{3!} \frac{\delta}{\delta J_1} \frac{\delta}{\delta J_2} \frac{\delta}{\delta J_3}\frac{\delta}{\delta J_4} \Big(D_{xz}J_{z} D_{xv}J_v D_{xw}J_w e^{-\frac{1}{2} J_z D_{zy}J_y} \partial^2 (D_{xy} J_y)\Big) \\
&= \frac{gh}{3!} \frac{\delta}{\delta J_1} \frac{\delta}{\delta J_2} \frac{\delta}{\delta J_3} \Big[3 D_{x4} D_{xv}J_v D_{xw} J_w \partial^2 (D_{xy} J_y) + D_{xz} J_z D_{xv} J_v D_{xw}J_w \partial^2 (D_{x4})\Big] e^{-\frac{1}{2}J_{z}D_{zy}J_{y}} \\
&= \frac{gh}{3!} \frac{\delta}{\delta J_1} \frac{\delta}{\delta J_2} \Big[3! D_{x4} D_{x3} D_{xw} J_w \partial^2(D_{xy} J_y) + 3 D_{x4} D_{xv} J_v D_{xw}J_w \partial^2 D_{x3} + 3 D_{x3}D_{xv} J_v D_{xw} J_w \partial^2 D_{x4}\Big]e^{-\frac{1}{2}J_z D_{zy}J_y} \\
&= \frac{gh}{3!} \frac{\delta}{\delta J_1} \Big[3! D_{x4} D_{x3} D_{x2} \partial^2 (D_{xy}J_y) + 3! D_{x4} D_{x3} D_{xw} J_w \partial^2 D_{x2} + 3! D_{x4} D_{x2} D_{xw} J_w \partial^2 D_{x3} + 3! D_{x3}D_{x2}D_{xw}J_w \partial^2 D_{x4}\Big]e^{-\frac{1}{2}J_z D_{zy}J_y} \\
&= gh [D_{x4} D_{x3} D_{x2} \partial^2 D_{x1} + D_{x4} D_{x3} D_{x1} \partial^2 D_{x2} + D_{x4}D_{x2}D_{x1}\partial^2 D_{x3} + D_{x3} D_{x2} D_{x1} \partial^2 D_{x4}]. \\
\end{align}
In the ordinary notation we find thus that:
\begin{equation}
\langle \Omega|T\{\phi(x_1)\phi(x_2)\phi(x_3)\phi(x_4)\}|\Omega \rangle = igh \int
d^4 x [D(x-x_4) D(x-x_3) D(x-x_2) \partial^2 D(x-x_1) + D(x-x_4) D(x-x_3) D(x-x_1) \partial^2 D(x-x_2) + D(x-x_4)D(x-x_2)D(x-x_1)\partial^2 D(x-x_3) + D(x-x_3) D(x-x_2) D(x-x_1) \partial^2 D(x-x_4)].
\end{equation}
Recall that the propagator is given by:
\begin{equation}
D(x-y) = \int \frac{d^4 k}{(2\pi)^4} \frac{i}{k^2 - m^2} e^{-ik \cdot x}. \\
\end{equation}
Therefore it is now easy to see that the vertex is given by:
\begin{equation}
V_{\phi \phi \phi \phi} = -igh [p_1^2 + p_2^2 + p_3^2 + p_4^2],
\end{equation}
where $p_1 + p_2 + p_3 + p_4 =0$ (so all momenta are incoming) since a delta function of the type $\delta(p_1 + p_2 + p_3 + p_4)$ will appear if you would compute the correlation function explicitly.
