In the first work of this series, we adopt an AT2017gfo-like viewing-angle-dependent kilonova model and the standard afterglow model with lightcurve distribution based on the cosmological short gammaray bursts afterglows to simulate the luminosity functions and color evolution of both kilonovae and optical afterglow emissions from binary neutron star mergers. Only a small fraction of (∼ 20%) on-axis afterglows are dimmer than the associated kilonovae at the kilonovae peak time. At a very large viewing angle, e.g., sin θ v 0.40, most kilonovae at the peak time would be much brighter than associated afterglows. Therefore, a large fractional of detectable kilonovae would be significantly polluted by associated afterglow emissions. We find that the maximum possible apparent magnitude for cosmological kilonovae is ∼ 27 − 28 mag. Afterglow emission has a wider luminosity function compared with kilonova emission. At brightness dimmer than ∼ 25 − 26 mag, according to their luminosity functions, the number of afterglows is much larger than that of kilonovae. Since the search depth of almost all the present and foreseeable survey projects is < 25 mag, the number of afterglow events detected via serendipitous observations would be much higher than that of kilonova events, consistent with the current observations. We find that it may be difficult to use the fading rate in a single band to directly identify kilonovae and afterglows among various fast-evolving transients by serendipitous surveys, especially if the observations have missed the peak time. However, the color evolution between optical and infrared bands can identify them, since their color evolution patterns are unique compared with those of other fast-evolving transients.