$mathit{Herschel}$ extragalactic surveys offer a unique opportunity to efficiently select a significant number of rare and massive dusty objects, and thus gain insight into the prodigious star-forming activity that takes place in the very distant Universe. To search for $zgeq4$ dusty star-forming galaxies, in this work we consider red SPIRE objects with fluxes rising from 250 $mu$m to $500:mu$m (so-called 500 $mu$m-risers). We aim to implement a novel method to obtain a statistical sample of 500 $mu$m-risers and fully evaluate our selection inspecting different models of galaxy evolution. We consider one of the largest and deepest ${it Herschel}$ surveys, the Herschel Virgo Cluster Survey. We develop a novel selection algorithm which links the source extraction and spectral energy distribution fitting. We select 133 500 $mu$m-risers over 55 deg$^{2}$, imposing the criteria: $S_{500}>S_{350}>S_{250}$, $S_{250}>13.2$ mJy and $S_{500}>$30 mJy. Differential number counts are in a fairly good agreement with models, displaying better match than other existing samples. In order to interpret the statistical properties of selected sources, which has been proven as a very challenging task due the complexity of observed artefacts, we make end-to-end simulations including physical clustering and lensing. The estimated fraction of strongly lensed sources is $24^{+6}_{-5}%$ based on models. We present the faintest known statistical sample of 500 $mu$m-risers and show that noise and strong lensing have crucial impact on measured counts and redshift distribution of selected sources. We estimate the flux-corrected star formation rate density at $4<z<5$ with the 500 $mu$m-risers and found it close to the total value measured in far-infrared. It indicates that colour selection is not a limiting effect to search for the most massive, dusty $z>4$ sources.