Consider the unit sphere, which can either be described by $x^2+y^2+z^2=1$ or by the equation $r(\theta,\phi)=1$, where $(r,\theta,\phi)$ are spherical polar coordinates.
I define a deformed sphere by the equation, $r(\theta,\phi) = 1+\delta r(\theta,\phi)$, where $\delta r (\theta,\phi)$ is a small smooth deformation of the surface. How would I go about writing a formal expression for the surface area of this deformed sphere, in a perturbative expansion in $\delta r$.
There are two conventions for spherical polar coordinates. For concreteness, I will use the physicist version here:
$$[0,2\pi] \times [0,2\pi) \ni (\theta,\phi) \quad\mapsto\quad \hat{n}(\theta,\phi) = (\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta) \in \mathbb{R}^3$$
When we deform the unit sphere slightly, we can continue to use $(\theta,\phi)$ as a parametrization of the deformed sphere. The position vector corresponding to $(\theta,\phi)$ will have the form
$$\vec{x}(\theta,\phi) = r(\theta,\phi) \hat{n}(\theta,\phi) = ( 1 + \epsilon \eta(\theta,\phi) )\hat{n}(\theta,\phi)$$ where $\epsilon$ is a small parameter controling the amount of deformation.
For any function $f$, let $f_\theta$ and $f_\phi$ be the partial derivative of $f$ with respect to $\theta$ and $\phi$ respectively. The area element of the deformed sphere is given by the formula
$$d\sigma = | \vec{x}_\theta \times \vec{x}_\phi | d\theta d\phi$$
Notice
$$\begin{cases} \vec{x}_\theta = r_\theta \hat{n} + r \hat{n}_\theta\\ \vec{x}_\phi = r_\phi \hat{n} + r \hat{n}_\phi\\ |\hat{n}| = |\hat{n}_\theta| = 1, |\hat{n}_\phi| = \sin\theta\\ \hat{n} \cdot \hat{n}_\theta = \hat{n} \cdot \hat{n}_\phi = \hat{n}_\theta \cdot \hat{n}_\phi = 0 \end{cases} \quad\implies\quad \begin{cases} |\vec{x}_\theta|^2 &= |r_\theta|^2 + r^2\\ |\vec{x}_\phi|^2 &= |r_\phi|^2 + r^2\sin^2\theta\\ \vec{x}_\theta \cdot \vec{x}_\phi &= r_\theta r_\phi \end{cases} $$ We have $$|\vec{x}_\theta \times \vec{x}_\phi| = \sqrt{ |\vec{x}_\theta|^2 |\vec{x}_\phi|^2 - |\vec{x}_\theta \cdot \vec{x}_\phi|^2} = r\sin\theta \sqrt{ r^2 + r_\theta^2 + \frac{r_\phi^2}{\sin^2\theta} } $$ This leads to $$d\sigma = \sin\theta d\theta d\phi \left[ (1 + \epsilon \eta )\sqrt{ (1 + \epsilon \eta)^2 + \epsilon^2 \Delta }\right] \quad\text{ where }\quad \Delta = \eta_\theta^2 + \frac{\eta_\phi^2}{\sin^2\theta} $$
Throw the RHS to a CAS and ask it to Taylor expand in $\epsilon$, we get
$$\frac{d\sigma}{\sin\theta d\theta d\phi} = 1 + 2\eta\epsilon + \left(\eta^2 + \frac{\Delta}{2}\right)\epsilon^2 - \frac{\Delta^2}{8} \epsilon^4 + \frac{\eta\Delta^2}{4} \epsilon^5 + \cdots $$ If you need more terms, you can ask a CAS to crank that out for you.