Если вам нужен именно вывод... Я бы вывел следующим образом:
Будем считать что функция плотности

задана на поверхности

и за ось примем

. Площадь сферы обозначим

Переходим к цилиндрическим координатам:

Мысленно разбиваем диаметр вдоль оси

на

равных отрезков и получаем соответственно разбиение сферы

на

сферических поясов (и шапочки на полюсах) одинаковой высоты

. Далее проводим

меридианов и соответственно разбиваем каждый сферический пояс на сферические "квадраты" ( в полюсах- треугольники) каждый с угловой мерой

. Получаем в общей сложности

кусочков.

В каждом кусочке /элементе сферы - выбираем точку с координатами

, и находим значение плотности

Момент каждого кусочка

относительно оси

приблизительно равен
Чтобы найти площадь

кривой поверхности, вообще говоря необходимо привлекать частные производные, но в нашем случае (сферы) спасает теорема, утверждающая, что если сферические пояса имеют одинаковую высоту, то их площади равны. (Доказывается средствами матана 2 курса или просто смотри в инете ). Поэтому равны и площади всех элементов нашего разбиения сферы

.
Далее все это суммируем по обоим индексам, устремляем

к бесконечности и получаем двойной интеграл по поверхности как советовали в самом начале.

Ну а если, как частный случай, масса распределена равномерно, то вытаскивая

за интегралы и умножив на

получаем

где

- масса сферы. Проинтегрировав, получите ответ данный
Архиповым.