Code and carrier phase based short baseline GPS positioning: computational aspects

A recursive least squares algorithm is presented for short baseline GPS positioning using both carrier phase and code measurements. We take advantage of the structure of the problem to make the algorithm computationally efficient and use orthogonal transformations to ensure that the algorithm is numerically reliable. Details are given for computing position estimates and error covariance matrices with possible satellite rising and setting. Real data test results suggest our algorithm is effective.This research was supported by NSERC of Canada Grant RGPIN217191–99, FCAR of Quebec Grant 2001-NC-66487, and NSERCGEOIDE Network Project ENV#14 for Xiao-Wen Chang, and by NSERC of Canada Grant RGPIN9236–01 for Christopher C. Paige.An erratum to this article can be found at     

A consistency test of airborne GPS using multiple monitor stations

In October 1990, several airborne GPS tests were conducted in the Ottawa region by the Canada Centre for Surveying (CCS) and the Canada Centre for Remote Sensing (CCRS). Ashtech XII receivers were located at up to three monitor stations with baseline lengths to the aircraft ranging from 1–200 km. Approximately two hours of airborne data, collected at a 2 Hz rate, were available for each of the three test days. Post-processing of the differential data was done using the University of Calgary's SEMIKIN package which utilizes a Kalman filter algorithm to estimate both the remote receiver's position and velocity. Comparisons were made between the aircraft position and velocity determined from each of the monitor stations to assess the consistency of differential GPS when different reference stations are used. Results show that the degree of consistency is dependent upon the distance to the monitor stations. Agreement at the decimetre-level is achieved in position when the baseline lengths are within 100 km. Agreement in velocity is usually better than 1 cm s^–1 (RMS).     

Excitation of non-atmospheric polar motion by the migration of the Pacific Warm Pool

Changes in the annual variation of the Earths polar motion are found to be largely caused by the variation of the atmospheric angular momentum (AAM). Recent simulation results of oceanic general circulation models further suggest global oceanic effects on the annual polar motion in addition to the atmosphere. In comparison with previous model studies of global oceanic effects, this research particularly singles out a large-scale ocean anomaly and investigates its effect on the annual polar motion, determined from satellite observations of the movement of the Western Pacific Warm Pool (WPWP). Although the scale of the warm pool is much smaller than that of the solid Earth, analysis of the non-atmospheric polar motion excitation has shown that the WPWP contributes non-negligibly to the annual polar motion. The analysis consists of over 30 years of WPWP data (1970–2000) and shows values of polar motion excitation for the x-component of (2.5 mas, –79°) and for the y-component of (0.6 mas, 173°). Comparison of this result with the total geodetic non-atmospheric polar motion excitation of (10.3 mas, 59°) for the x-component and (10.6 mas, 62°) for the y-component shows the significance of the WPWP. Changes in the Earths polar motion have attracted significant attention, not only because it is an important geodetic issue, but also because it has significant value as a global measure of variations within the hydrosphere, atmosphere, cryosphere, and solid Earth, and hence global changes.Tel: 86–21–64386191 Fax: 86–21–64384618Acknowledgments. The authors are grateful to Dr. R. Gross (JPL) and two anonymous reviewers for providing invaluable comments. They also thank Dr. J.L. Chen (CSR) for helpful discussions. Y. Zhou, D. Zheng and X. Liao were supported by the National Natural Science Foundation of China (10273018, 10133010) and Key Project of Chinese Academy of Sciences (KJCX2-SW-T1). X-H. Yan was supported by the National Aeronautics and Space Administration (NASA) through Grant NAG5–12745, and by the National Science Foundation (NSF) through the Presidential Faculty Fellow award to X-H. Yan (OCE-9453499). W.T. Liu was supported by the NASA Physical Oceanography Program.     

Wavelet evaluation of the Stokes and Vening Meinesz integrals

The wavelet transform is a powerful tool in evaluating some singular geodetic integrals. Due to its localization properties in both the time (space) and frequency (scale) domains, and because the kernels of some geodetic integrals have singular points and decay smoothly and quickly away from the singularities, many wavelet transform coefficients of the kernels become zeros or negligible, and only a small number of wavelet transform coefficients are significant. It is thus possible to significantly compress the kernels of these integrals on a wavelet basis by neglecting the zero coefficients and the small coefficients below a certain threshold. Therefore, wavelets provide a convenient way of efficiently evaluating these integrals in terms of fast computation and savings of computer memory. A modified algorithm for the wavelet evaluation of Stokes' integral is presented. The same modified algorithm is applied to the evaluation of the Vening Meinesz integral, whose kernel has a stronger singularity than does Stokes' kernel. Numerical examples illustrate the efficiency and accuracy of the wavelet methods.Acknowledgments.The author express their sincere thanks to Dr. Salamonwicz for providing his PhD thesis. E-mail correspondence between the authors and Dr. Barthelmes and Dr. Benciolini contributed to the work. R. Benciolini and the other two anonymous reviewers are thanked for their constructive comments. Support for this research was provided by research grants to Dr. Sideris from the Natural Sciences and Engineering Reserch Council of Canada (NSERC) and the Geomatics for Informed Decisions (GEOIDE) Network of Centres of Excellence. The MATLAB Wavelet Toolbox package was used as the platform to develop the software in this project.     

Applications of downward-continuation in gravimetric geoid modeling: case studies in Western Canada

The objective of this study is to evaluate two approaches, which use different representations of the Earth’s gravity field for downward continuation (DC), for determining Helmert gravity anomalies on the geoid. The accuracy of these anomalies is validated by 1) analyzing conformity of the two approaches; and 2) converting them to geoid heights and comparing the resulting values to GPS-leveling data. The first approach (A) consists of evaluating Helmert anomalies at the topography and downward-continuing them to the geoid. The second approach (B) downward-continues refined Bouguer anomalies to the geoid and transforms them to Helmert anomalies by adding the condensed topographical effect. Approach A is sensitive to the DC because of the roughness of the Helmert gravity field. The DC effect on the geoid can reach up to 2 m in Western Canada when the Stokes kernel is used to convert gravity anomalies to geoid heights. Furthermore, Poisson’s equation for DC provides better numerical results than Moritz’s equation when the resulting geoid models are validated against the GPS-leveling. On the contrary, approach B is significantly less sensitive to the DC because of the smoothness of the refined Bouguer gravity field. In this case, the DC (Poisson’s and Moritz’s) contributes only at the decimeter level to the geoid model in Western Canada. The maximum difference between the geoid models from approaches A and B is about 5 cm in the region of interest. The differences may result from errors in the DC such as numerical instability. The standard deviations of the h−H−N for both approaches are about 8 cm at the 664 GPS-leveling validation stations in Western Canada.     

Economic convergence or divergence? Modeling the interregional dynamics of EU regions, 1985–1999

In this paper we start from a continuous time framework derived from the classical predator-prey model in order to analyze the recent dynamics of regional evolution in the EU. The model describes a system of interrelated units obeying a complex functional dynamics that at any moment may encompass divergent forces. After briefly reviewing the modeling framework presented elsewhere (Arbia and Paelinck 2002), we consider the dynamics of per-capita income in 119 NUTS-2 European regions in the years 1985–1999, and fit the model using Simultaneous Dynamic Least Squares (Paelinck 1996b). The model is shown to fit the empirical data well. The – potential divergence – leads up to a formidable coordination task for public authorities responsible for regional policies, as they are supposed – if their aim is absolute convergence in an economic sense, meaning equalizing incomes per head – to produce at any time a consistent set of regional policies for all the regions concerned, and not only for the less developed ones.     

Estimating forest biomass using small footprint LiDAR data: An individual tree-based approach that incorporates training data

A new individual tree-based algorithm for determining forest biomass using small footprint LiDAR data was developed and tested. This algorithm combines computer vision and optimization techniques to become the first training data-based algorithm specifically designed for processing forest LiDAR data. The computer vision portion of the algorithm uses generic properties of trees in small footprint LiDAR canopy height models (CHMs) to locate trees and find their crown boundaries and heights. The ways in which these generic properties are used for a specific scene and image type is dependent on 11 parameters, nine of which are set using training data and the Nelder–Mead simplex optimization procedure. Training data consist of small sections of the LiDAR data and corresponding ground data. After training, the biomass present in areas without ground measurements is determined by developing a regression equation between properties derived from the LiDAR data of the training stands and biomass, and then applying the equation to the new areas. A first test of this technique was performed using 25 plots (radius = 15 m) in a loblolly pine plantation in central Virginia, USA (37.42N, 78.68W) that was not intensively managed, together with corresponding data from a LiDAR canopy height model (resolution = 0.5 m). Results show correlations (r) between actual and predicted aboveground biomass ranging between 0.59 and 0.82, and RMSEs between 13.6 and 140.4 t/ha depending on the selection of training and testing plots, and the minimum diameter at breast height (7 or 10 cm) of trees included in the biomass estimate. Correlations between LiDAR-derived plot density estimates were low (0.22 ≤ r ≤ 0.56) but generally significant (at a 95% confidence level in most cases, based on a one tailed test), suggesting that the program is able to properly identify trees. Based on the results it is concluded that the validation of the first training data-based algorithm for determining forest biomass using small footprint LiDAR data was a success, and future refinement and testing are merited.     

Radiometric correction of hemispherical images

Light scattering in fore-optics of imaging radiometric instruments degrades images, the signal of dark pixels is enhanced and that of bright pixels is reduced. A procedure of restoring hemispherical images obtained with a hemispheric view CCD-radiometer is suggested. The circular hemispheric image is projected onto the sphere and Wiener filtering on the sphere is performed. For building the Wiener filter the point spread function of the instrument was measured in the laboratory, and test images of (partially) controlled radiance were measured. The procedure is applied for the correction of hemispherical radiometric images taken under forest canopy. In non-corrected near infrared images the signal of sky pixels in forest canopy gaps near zenith is enhanced 1.4–2.0 times. The applied correction reduces this error to the level of 1.1–1.2. In the red channel these figures are 0.8–0.9 and 0.97–0.99, respectively.     

Testing satellite and ground thermal imaging of low-temperature fumarolic fields: The dormant Nisyros Volcano (Greece)

The Nisyros Volcano (Greece) was monitored by satellite and ground thermal imaging during the period 2000–2002. Three night-scheduled Landsat-7 ETM⁺ thermal (band 6) images of Nisyros Island were processed to obtain land surface temperature. Ground temperature data were also collected during one of the satellite overpasses. Processed results involving orthorectification and 3-D atmospheric correction clearly show the existence of a thermal anomaly inside the Nisyros Caldera. This anomaly is associated mainly with the largest hydrothermal craters and has land surface temperatures 5–10 °C warmer than its surroundings. The ground temperature generally increased by about 4 °C inside the main crater over the period 2000–2002. Ground thermal images of the hydrothermal Stephanos Crater were also collected in 2002 using a portable thermal infrared camera. These images were calibrated to ground temperature data and orthorectified. A difference of about 0–2 °C was observed between the ground thermal images and the ground temperature data. The overall study demonstrates that satellite remote sensing of low-temperature fumarolic fields within calderas can provide a reliable long-term monitoring tool of dormant volcanoes that have the potential to reactivate. Similarly, a portable thermo-imager can easily be deployed for real-time monitoring using telemetric data transfer. The operational costs for both systems are relatively low for an early warning system.     

Comparison of resection–intersection algorithms and projection geometries in radiostereometry

Three resection–intersection algorithms were applied to simulated projections and clinical data from radiostereometric patients. On simulated data, the more advanced bundle-adjustment-based algorithms outperformed the classical Selvik algorithm, even if the error reductions were small for some parameters. On clinical data, the results were inconclusive.The two different projection geometries had a much larger influence on the error size and distribution. For the biplanar configuration, the position and motion errors were small and almost isotropic. For the uniplanar configuration, the position errors were comparably high and anisotropic, but still resulted in a high accuracy for some motion parameters at the expense of others.The simplified resection–intersection algorithm by Selvik may still be considered a good and robust algorithm for radiostereometry. More studies will have to be performed to find out how the theoretical advantages of the bundle methods can be utilized in clinical radiostereometry.     

The ADS40 Vaihingen/Enz geometric performance test

This paper presents the main results of a comprehensive ADS40 performance analysis conducted at the Vaihingen/Enz test field. As such it represents one example of an independent in-flight performance study for one of the new and commercially operational digital airborne camera systems. Based on a large number of well coordinated and defined object points, which served as independent check points, the absolute geometric accuracy of the ADS40 from true operational data has been verified. Empirical analysis of data from flying heights ranging from 1500 m to 4000 m proved the ADS40 geometric accuracy to be in the range of 1–2 μm at image scale for horizontal coordinates and 0.03–0.05‰ of the flying height for vertical components. This is fully within specification for airborne imaging.     

Scale factor calibration of a superconducting gravimeter at Esashi Station, Japan, using absolute gravity measurements

The scale factor of a superconducting gravimeter (SG) at the Esashi Earth Tides Station, Japan, was revised by repeating co-located absolute gravity measurements with an FG5 gravimeter. Although the calibration results from the absolute gravimeter (AG) show an apparent secular change in the scale factor of the SG (0.4% for the period 1993–2002), the relative scale factors, which are determined by tidal analysis with the response method, indicate that it has changed by no more than 0.01% during the above period. If the mean scale factor over the 10 years is adopted, a value of –56.082±0.029 Gal/V (1 Gal =10^–8 m s^–2) is obtained, which is about 0.4% smaller than that used in the global geodynamics project (GGP) database. Based on this newly determined scale factor, the tidal gravity factors at Esashi have been re-estimated. The observed tidal factors, corrected for the ocean tide effects with recent models, indicate that the theoretical gravity factors for an inelastic Earth model are more consistent with the observations than are those for an elastic model.     

Effect of the atmospheric pressure on surface displacements

Summary Atmospheric pressure variations with periods of some days and months can be considered as loading functions on the Earth's surface and can induce quasi-periodic surface deformations. The influence of such surface displacements is calculated by performing a convolution sum between the mass loading Green's functions and the local and regional barometric pressure data (geographically distribution in a 1° × 1° grid system extending to more than 1000km). The results for 5 stations in Europe show that the average values reach about 22.9–30.2mm depending on the ocean response: the inverted or non-inverted barometer ocean model. The corresponding admittances are 0.576–0.758mm/mbar respectively. The horizontal displacements are not negligible but always smaller. The magnitudes are about 2–3mm for East-West component and 0.5–1.0mm for North-South component.The results of the dependence on the lateral extension of the pressure load show that the admittance for radial displacement varies from 0.250mm/mbar for a column load of 100km radius to 0.539mm/mbar for a column load of more than 1000km extension. It means that the main contribution of the loads comes from the horizontal distribution of the air pressure in a broad region.The time dependent effects of the atmospheric pressure are also computed with the two-coefficient correction equations provided by Rabbel & Zschau (1985) using ground pressure data in a 1.125° × 1.125° grid system. The computations demonstrate that the results are in good agreement with those obtained with a convolution sum. It shows that this method can provide us with a good approximation for vertical displacement caused by the deformation of the Earth.     

Assessing the spatial distribution of crop areas using a cross-entropy method

While crop production statistics are reported on a geopolitical – often national – basis, we often need to know, for example, the status of production or productivity within specific sub-regions, watersheds, or agro-ecological zones. Such re-aggregations are typically made using expert judgments or simple area-weighting rules. We describe a new, entropy-based approach to the plausible estimates of the spatial distribution of crop areas. Using this approach tabular crop production statistics are blended judiciously with an array of other secondary data to assess the areas of specific crops within individual ‘pixels’—typically 25–100 km² in size. The information utilized includes crop production statistics, farming system characterization, satellite-based interpretation of land cover, biophysical crop suitability assessments, and population density. An application is presented in which Brazilian state level production statistics are used to generate pixel level crop area data for eight crops. To validate the spatial allocation we aggregated the pixel estimates to obtain synthetic estimates of municipality level areas in Brazil, and compared those estimates with actual municipality statistics. The approach produced extremely promising results. We then examined the robustness of these results compared to simplified approaches to spatializing crop production statistics and showed that, while computationally intensive, the cross-entropy method does provide more reliable spatial allocations.     

Evaluation of approaches for AWiFS multi-date registration

AWiFS onboard IRS-P6 belongs to the category of high-repetivity sensors based on large swath, but with ground trace based on narrow-swath sensor (LISS-III). This is useful for cloud removal as well as vegetation phenology studies. Such multi-date analysis has a prerequisite of accurate multi-date registration. This study investigates the accuracy of multi-date registration over a mixed plain and hilly terrain in northern India (29–31°N latitude and 77.5–79.5°E longitude; 200–4000 m.a.s.l.). Simple polynomial rectification, multi-date registration using ortho-correction technique on standard product (level-2) and radiometric product (level-1) as a function of number of ground control points (GCPs) and external Digital Elevation Model (DEM) were investigated. The results indicated that ortho-rectification on level-1 product provided better accuracy in comparison to simple rectification and ortho-rectification on level-2 product.     

Gravimetric examination of Hagia Sophia's subsurface structure

The subsurface structure of Hagia Sophia, one of the oldest sacred monuments in the world built between 532–537 under the reign of Justinian in today's Istanbul, has been investigated by using two relative LaCoste-Romberg gravimeters in order to detect hidden cavities which have also served as earthquake dampers in similar constructions. On the building's ground floor a grid of 100 points with a grid size of about 4.m was measured. The mean gravimetric point error was ± 3.10^–8 ms^–2. The result of the examination is that cavities were not detected in the inner central part of Hagia Sophia with a larger diameter than 8.m down to a depth of about 20.m, and Hagia Sophia's foundation was found to be a slope of natural rock with a downward inclination to the East that has a small crest symmetrical to the building's East-West axis.     

A bilinear approximation of the surface relief in terrain correction computations

A new method is presented for the computation of the gravitational attraction of topographic masses when their height information is given on a regular grid. It is shown that the representation of the terrain relief by means of a bilinear surface not only offers a serious alternative to the polyhedra modeling, but also approaches even more smoothly the continuous reality. Inserting a bilinear approximation into the known scheme of deriving closed analytical expressions for the potential and its first-order derivatives for an arbitrarily shaped polyhedron leads to a one-dimensional integration with – apparently – no analytical solution. However, due to the high degree of smoothness of the integrand function, the numerical computation of this integral is very efficient. Numerical tests using synthetic data and a densely sampled digital terrain model in the Bavarian Alps prove that the new method is comparable to or even faster than a terrain modeling using polyhedra.     

Ellipsoidal terrain correction based on multi-cylindrical equal-area map projection of the reference ellipsoid

An operational algorithm for computation of terrain correction (or local gravity field modeling) based on application of closed-form solution of the Newton integral in terms of Cartesian coordinates in multi-cylindrical equal-area map projection of the reference ellipsoid is presented. Multi-cylindrical equal-area map projection of the reference ellipsoid has been derived and is described in detail for the first time. Ellipsoidal mass elements with various sizes on the surface of the reference ellipsoid are selected and the gravitational potential and vector of gravitational intensity (i.e. gravitational acceleration) of the mass elements are computed via numerical solution of the Newton integral in terms of geodetic coordinates {,,h}. Four base- edge points of the ellipsoidal mass elements are transformed into a multi-cylindrical equal-area map projection surface to build Cartesian mass elements by associating the height of the corresponding ellipsoidal mass elements to the transformed area elements. Using the closed-form solution of the Newton integral in terms of Cartesian coordinates, the gravitational potential and vector of gravitational intensity of the transformed Cartesian mass elements are computed and compared with those of the numerical solution of the Newton integral for the ellipsoidal mass elements in terms of geodetic coordinates. Numerical tests indicate that the difference between the two computations, i.e. numerical solution of the Newton integral for ellipsoidal mass elements in terms of geodetic coordinates and closed-form solution of the Newton integral in terms of Cartesian coordinates, in a multi-cylindrical equal-area map projection, is less than 1.6×10^–8 m²/s² for a mass element with a cross section area of 10×10 m and a height of 10,000 m. For a mass element with a cross section area of 1×1 km and a height of 10,000 m the difference is less than 1.5×10^–4m²/s². Since 1.5× 10^–4 m²/s² is equivalent to 1.5×10^–5m in the vertical direction, it can be concluded that a method for terrain correction (or local gravity field modeling) based on closed-form solution of the Newton integral in terms of Cartesian coordinates of a multi-cylindrical equal-area map projection of the reference ellipsoid has been developed which has the accuracy of terrain correction (or local gravity field modeling) based on the Newton integral in terms of ellipsoidal coordinates.Acknowledgments. This research has been financially supported by the University of Tehran based on grant number 621/4/859. This support is gratefully acknowledged. The authors are also grateful for the comments and corrections made to the initial version of the paper by Dr. S. Petrovic from GFZ Potsdam and the other two anonymous reviewers. Their comments helped to improve the structure of the paper significantly.     

Correction of laser scanning intensity data: Data and model-driven approaches

Most airborne and terrestrial laser scanning systems additionally record the received signal intensity for each measurement. Multiple studies show the potential of this intensity value for a great variety of applications (e.g. strip adjustment, forestry, glaciology), but also state problems if using the original recorded values. Three main factors, a) spherical loss, b) topographic and c) atmospheric effects, influence the backscatter of the emitted laser power, which leads to a noticeably heterogeneous representation of the received power. This paper describes two different methods for correcting the laser scanning intensity data for these known influences resulting in a value proportional to the reflectance of the scanned surface. The first approach – data-driven correction – uses predefined homogeneous areas to empirically estimate the best parameters (least-squares adjustment) for a given global correction function accounting for all range-dependent influences. The second approach – model-driven correction – corrects each intensity independently based on the physical principle of radar systems. The evaluation of both methods, based on homogeneous reflecting areas acquired at different heights in different missions, indicates a clear reduction of intensity variation, to 1/3.5 of the original variation, and offsets between flight strips to 1/10. The presented correction methods establish a great potential for laser scanning intensity to be used for surface classification and multi-temporal analyses.