We report results on the studies of the $e^+e^-\to B_s^{(*)}\bar{B}_s^{(*)}$
processes. The results are based on a $121.4$ fb$^{-1}$ data sample collected
with the Belle detector at the center-of-mass energy near the $\Upsilon(10860)$
peak and $16.4$ fb$^{-1}$ of data collected at 19 energy points in the range
from 10.77 to 11.02 GeV. We observe a clear $e^+e^-\to\Upsilon(10860)\to
B_s^{(*)}\bar{B}_s^{(*)}$ signal, with no statistically significant signal of
$e^+e^-\to \Upsilon(11020)\to B_s^{(*)}\bar{B}_s^{(*)}$. The relative
production ratio of $B_s^*\bar{B}_s^*$, $B_s\bar{B}_s^{*}$, and $B_s\bar{B}_s$
final states at $\sqrt{s}=10.866$ GeV is measured to be $7:$
$0.856\pm0.106(stat.)\pm0.053(syst.):$
$0.645\pm0.094(stat.)^{+0.030}_{-0.033}(syst.)$. An angular analysis of the
$B_s^*\bar{B}_s^*$ final state produced at the $\Upsilon(10860)$ peak is also
performed.