I am thinking about a creative use of the Euler Beta function.
For now, just a symmetry trick related to the formula $\tan(\pi/4-u)=\frac{1-\tan u}{1+\tan u}$:
$$\begin{eqnarray*}I&=&\int_{0}^{1}\frac{\log(1+x)}{1+x^2}\,dx = \int_{0}^{\pi/4}\log(1+\tan\theta)\,d\theta\\&=&\int_{0}^{\pi/4}\log\left(1+\frac{1-\tan\theta}{1+\tan\theta}\right)\,d\theta=\frac{\pi}{4}\log 2-I\end{eqnarray*}$$
from which:
$$ I = \color{red}{\frac{\pi}{8}\log 2}.$$
I would recommend you to work through this paper which introduces the Gamma aswell as the Beta Function quite nicely and also provides a proof of their infamous relation.
However, the key theorem we are using will working with substitutions within double (or triple, etc.) integrals is a general substitution formula given by
$$\int_{\varphi(U)}f(\mathbf v)\mathrm d\mathbf v~=~\int_Uf(\varphi(\mathbf u))|\det(D\varphi)(\mathbf u)|\mathrm d\mathbf u\tag1$$
Here $\mathbf v$ and $\mathbf u$ are same dimensional vectors used to shorten things up and $\varphi$ is our substitution. What is really intersting about this formula is the differential on the RHS. $|\det(D\varphi)(\mathbf u)|$ refers to the absolute value of the so-called Jacobian Determinant, which is a matrix consisting of partial derivatives of every $\varphi(u_i)$ w.r.t. every $u_i$. Note that for the case that we are dealing with $1$-dimensional functions $(1)$ is just our well-known Integration By Substitution.
Let us get back to the given case. We want to compute the following integral using the subsitutions $x=vt$ and $y=v(1-t)$ (note that there is a little mistake within your post as pointed out by kimchi lover)
$$\int_0^\infty\int_0^\infty x^{m-1}y^{n-1}e^{-(x+y)}\mathrm dx\mathrm dy$$
What we need to compute now are $1.$ the absolute value of the Jacobian Determinant and $2.$ the new borders of integration. The first one is quite simple since we get that
$$|\det(D\varphi(v,t))|=\begin{vmatrix}\frac{\partial}{\partial v}(vt)&\frac{\partial}{\partial t}(vt)\\\frac{\partial}{\partial v}(v(1-t))&\frac{\partial}{\partial t}(v(1-t))\end{vmatrix}=\begin{vmatrix}t&v\\1-t&-v\end{vmatrix}=|-vt-(1-t)v|=|-v|=v$$
Thus, we know that $\mathrm dx\mathrm dy=v~\mathrm dv\mathrm dt$. Since $0<x<\infty$ and $0<y<\infty$ we get by the relation $v=x+y$ and $t=\frac x{x+y}$ that $0<v<\infty$ and $0<t<1$ as our new borders (to be honest this part is tricky sometimes).
Combining what we have accomplished before we finally obtain by $(1)$ that
$$\therefore~\int_0^\infty\int_0^\infty x^{m-1}y^{n-1}e^{-(x+y)}\mathrm dx\mathrm dy=\int_0^1\int_0^\infty (vt)^{m-1}(v(1-t))^{n-1}e^{-v}v~\mathrm dv\mathrm dt$$
From hereon it is rather easy to be finished, i.e. splitting up the double integral, applying the defintions of the Gamma and the Beta Function and you are done.
Best Answer
Deriving the Surface Area of a Sphere
Start with the $1$ dimensional integral $$ \int_{-\infty}^\infty e^{-\pi x^2}\,\mathrm{d}x=1\tag1 $$ and take the product of $n$ copies: $$ \int_{-\infty}^\infty\int_{-\infty}^\infty\cdots\int_{-\infty}^\infty e^{-\pi \left(x_1^2+x_2^2+\cdots+x_n^2\right)}\,\mathrm{d}x_1\,\mathrm{d}x_2\cdots\,\mathrm{d}x_n=1\tag2 $$ Since the integrand is constant over spherical shells, we can easily write it as a spherical integral $$ \int_0^\infty e^{-\pi r^2}\omega_{n-1}r^{n-1}\,\mathrm{d}r=1\tag3 $$ where $\omega_{n-1}$ is the surface area of the $n-1$ dimensional unit sphere.
Thus, we get $$ \begin{align} 1 &=\int_0^\infty e^{-\pi r^2}\omega_{n-1}r^{n-1}\,\mathrm{d}r\tag{4a}\\ &=\frac{\omega_{n-1}}2\int_0^\infty e^{-\pi r^2}r^{n-2}\,\mathrm{d}r^2\tag{4b}\\ &=\frac{\omega_{n-1}}2\int_0^\infty e^{-\pi r}r^{n/2-1}\,\mathrm{d}r\tag{4c}\\ &=\frac{\omega_{n-1}}{2\pi^{n/2}}\int_0^\infty e^{-r}r^{n/2-1}\,\mathrm{d}r\tag{4d}\\ &=\frac{\omega_{n-1}}{2\pi^{n/2}}\Gamma(n/2)\tag{4e} \end{align} $$ Explanation:
$\text{(4a)}$: copy $(3)$
$\text{(4b)}$: $r\,\mathrm{d}r=\frac12\mathrm{d}r^2$ and pull the constant out front
$\text{(4c)}$: substitute $r\mapsto r^{1/2}$
$\text{(4d)}$: substitute $r\mapsto r/\pi$
$\text{(4e)}$: definition of the Gamma function
Therefore, we get $$ \omega_{n-1}=\frac{2\pi^{n/2}}{\Gamma(n/2)}\tag5 $$
Application to the Question
We can write the integral from the question as an integral on spherical shells: $$ \begin{align} &\int_0^\infty\overbrace{(a+r)^me^{-\lambda r^{2k}}}^{\substack{\text{integrand on a}\\\text{shell of radius $r$}}}\overbrace{\omega_{n-1}r^{n-1\vphantom{r^2}}}^{\substack{\text{surface area}\vphantom{g}\\\text{of the shell}}}\,\mathrm{d}r\tag{6a}\\ &=\frac{\omega_{n-1}}{2k}\sum_{j=0}^m\binom{m}{j}\int_0^\infty a^{m-j}r^{\frac{j+n}{2k}-1}e^{-\lambda r}\,\mathrm{d}r\tag{6b}\\ &=\frac{\omega_{n-1}}{2k}\sum_{j=0}^m\binom{m}{j}a^{m-j}\lambda^{-\frac{j+n}{2k}}\int_0^\infty r^{\frac{j+n}{2k}-1}e^{-r}\,\mathrm{d}r\tag{6c}\\ &=\frac{\omega_{n-1}}{2k}\sum_{j=0}^m\binom{m}{j}a^{m-j}\lambda^{-\frac{j+n}{2k}}\Gamma\!\left(\frac{j+n}{2k}\right)\tag{6d} \end{align} $$ Explanation:
$\text{(6a)}$: integral from the question
$\text{(6b)}$: substitute $r\mapsto r^{\frac1{2k}}$ and apply the Binomial Theorem
$\text{(6c)}$: substitute $r\mapsto r/\lambda$
$\text{(6d)}$: definition of the Gamma function