diff --git a/skyfield/starlib.py b/skyfield/starlib.py
index 1090e1fc..6aec090b 100644
--- a/skyfield/starlib.py
+++ b/skyfield/starlib.py
@@ -117,13 +117,13 @@ def _observe_from_bcrs(self, observer):
         if len(position.shape) > 1:
             if len(observer.xyz.au.shape) > 1:
                 vector = position - observer.xyz.au
-                vel = (observer.velocity.au_per_d.T - velocity).T
+                vel = (velocity - observer.velocity.au_per_d.T).T
             else:
                 vector = (position.T - observer.xyz.au).T
-                vel = (observer.velocity.au_per_d - velocity.T).T
+                vel = (velocity.T - observer.velocity.au_per_d).T
         else:
             vector = position - observer.xyz.au
-            vel = (observer.velocity.au_per_d.T - velocity).T
+            vel = (velocity - observer.velocity.au_per_d.T).T
         distance = length_of(vector)
         light_time = distance / C_AUDAY
         if len(position.shape) > 1 and not t.shape:
diff --git a/skyfield/tests/test_stars.py b/skyfield/tests/test_stars.py
index a77e9533..3921a418 100644
--- a/skyfield/tests/test_stars.py
+++ b/skyfield/tests/test_stars.py
@@ -1,5 +1,17 @@
 from skyfield import api
 from skyfield.data.hipparcos import load_dataframe
+from skyfield.framelib import ICRS
+
+def test_astrometric_velocity():
+    ts = api.load.timescale()
+    for zero, jd in ((0.0, api.T0),
+                     (0.0, [api.T0, api.T0 + 1.0]),
+                     ([0.0, 0.0], api.T0)):
+        star = api.Star(ra_hours=zero, dec_degrees=zero,
+                        ra_mas_per_year=1000.0, parallax_mas=1000.0)
+        a = api.SSB.at(ts.tdb_jd(jd)).observe(star)
+        ra_rate = a.frame_latlon_and_rates(ICRS)[4]
+        assert (abs(ra_rate.mas.per_day - 1000.0 / 365.25) < 1e-9).all()
 
 def test_dataframe():
     with api.load.open('hip_main.dat.gz') as f: