Very precise kinematic or dynamic orbits based on measurements of the Global Navigation Satellite Systems (GNSS) are required to study sea level change and ice cover variations based on the observations of altimetry satellites, atmospheric sounding by GNSS occultation measurements or the detection of mass transports and the mass distribution in the Earth system by a precise determination of the stationary and time variable gravity field. The continuous and precise observation of the orbits of low flying satellites such as CHAMP and GRACE by the GNSS enabled the development of new gravity field determination techniques. The classical approach of satellite geodesy was based on the analysis of accumulated orbit perturbations of artificial satellites with different altitudes and orbit inclinations. This so-called differential orbit improvement technique required the analysis of rather long arcs of days to weeks; it was the adequate technique for satellite arcs poorly covered with observations, mainly precise laser ranging to satellites. Now a very dense coverage with observations of the low flying satellites is available independent from Earth based observation stations and there is no need to use very long arcs with its intrinsic problems. The new alternative gravity field recovery concepts developed in the last couple of years require precise kinematical or reduced dynamical orbits derived from the code and phase measurements. In this paper a new approach for an integrated kinematic-dynamic orbit determination of low flying satellites based on GNSS observations is presented. The short arcs of the low flying satellites are represented by a linear approximation function where the model parameters are also functions of the force function acting on the satellites. This allows the determination of orbits with different kinematic and dynamic orbit characteristics.