Two strongly interacting bosons are trapped in an incommensurate model, which is described as $\hat H = - J\sum\limits_j{} {\left( {\hat c_j^\dagger {{\hat c}_{j + 1}} + {\rm{h}}{\rm{.c}}{\rm{.}}} \right)} + 2\lambda \sum\limits_j{} {\dfrac{{\cos \left( {2{\text{π}}\alpha j} \right)}}{{1 - b\cos \left( {2{\text{π}}\alpha j} \right)}}} {\hat n_j} + \dfrac{U}{2}\sum\limits_j{} {{{\hat n}_j}\left( {{{\hat n}_j} - 1} \right)} ,$ where there exists no interaction, the system displays mobility edges at $b\varepsilon = 2(J - \lambda )$, which separates the extended regime from the localized one and b = 0 is the standard Aubry-André model. By applying the perturbation method to the third order in a strong interaction case, we can induce an effective Hamiltonian for bosonic pairs. In the small b case, the bosonic pairs present the mobility edges in a simple closed expression form $b\left( {\dfrac{{{E^2}}}{U} - E - \dfrac{4}{E}} \right) = - 4\left(\dfrac{1}{E} + \lambda \right)$, which is the central result of the paper. In order to identify our results numerically, we define a normalized participation ratio (NPR) $\eta (E)$ to discriminate between the extended properties of the many-body eigenvectors and the localized ones. In the thermodynamic limit, the NPR tends to 0 for a localized state, while it is finite for an extended state. The numerical calculations finely coincide with the analytic results for b = 0 and small b cases. Especially, for the b = 0 case, the mobility edges of the bosonic pairs are described as $\lambda = - 1/E$<