The saturation magnetisation is the maximal magnetisation that your system can exhibit. On simple lattice systems with not too complex interactions this is just the point where all your spins are aligned in some particular direction.
In order to link it to temperature and magnetic moments, there are several models available. I categorise them in two categories: mean field theories and fully quantum mechanical approach. In your case I think that a mean field approach can be useful to obtain the expression that you are after. I will sketch a way to obtain it now, and if needed will upload a more complete answer tonight.
Essentially you consider the statistical mechanics of magnetic moments in an external field H. The magnetisation will be given by the derivative of the free energy wrt your field (up to a minus sign to make the thermodynamics correct). Once you have this expression, you replace your field H with H + cM where c is some constant and M the magnetisation. cM is essentially a molecular field term introduced by Weiss: the field generated by the static atoms on the lattice is proportional to the magnetisation. You then obtain a transcendental equation which you can solve numerically. Within the correct limits you can Taylor expand one side to simplify it a bit further. Hope this helps.