A. M. Loghin, R.S Ajin, V. E. Oniga [615884]

A. M. Loghin, R.S Ajin, V. E. Oniga
The threeAdimensional geodetic networks adjustement automation using Matlab

15
THE THREE-DIMENSIONAL GEODETIC NETWORKS
ADJUSTEMENT AUTOMATION USING MATLAB

AnaVMaria LOGHIN –Ph.D Student: [anonimizat]., „Gheorghe Asa chi” Technical University of Iasi,
Faculty of Hydrotechnical Engineering, Geodesy and Environmental Engineering,
[anonimizat]
Ajin.R.S – Independent Researcher, ajinares@gmail.c om
Valeria Ersilia ONIGA –Lecturer, Ph.D Eng., „Gheorg he Asachi” Technical University of
Iasi, Faculty of Hydrotechnical Engineering, Geodes y and Environmental Engineering,
Department of Terrestrial Measurements and Cadastre , [anonimizat]

Abstract: A threeVdimensional geodetic network is composed by a system of points
located on earth’s surface, on an engineering work territory, points which have a known
position in an unitary reference coordinate system. Geodetic networks play a significant role
in the attainable accuracy of surveying and mapping applications. This paper presents a
threeVdimensional geodetic network in form of a clo sedVloop traverse, which begins and ends
at the same point station. In order to determine th e traverse station points coordinates, the
least squares method of indirect observations was u sed. The primary objective of this paper is
writing a Matlab program with procedures and functi ons, which automatically calculates the
unknown points coordinates, respectively their accu racy, thus streamlining the measurements
processing.
Keywords: threeVdimensional geodetic network, closedVloop tra verse, least squares
method, Matlab

1. Introduction

It is already known that no observation made to a t hing or a phenomena will never be
exact. Every observation contains errors. These sta tements are fundamental and accepted
universally. Therefore, surveyors, who are measurem ent specialists, should have a thorough
understanding of errors. They must be familiar with the different types of errors, their sources
and their expected magnitudes and they must be able to adopt procedures for reducing error
sizes when making their measurements and to analyze and adjust their data.
In recent years, the least square method of adjusti ng spatial data has been rapidly
gaining popularity as the method used for analyzing and adjusting surveying data. This fact
should not be surprising, because this method is th e most rigourous adjustment procedure
available [1].
Geodetic surveying is concerned with two different types of surveying: the first one
for heights and the second one for horizontal eleme nts. In recent years, with the technology
evolution, there was developed and implemented the concept of threeAdimensional geodetic
network, that is a fusion between horizontal and ve rtical parts. Therefore, a threeAdimensional
geodetic network includes both types of measurement s: for heights and also for the horizontal
elements being composed by a system of points, loca ted on the earth’s surface, on an
engineering work territory, points which have a kno wn position in an unitary reference
coordinate system.

„1 Decembrie 1918” University of Alba Iulia RevCA D 19/2015

16
With the development of the technology of Electroni c Distance measurement,
traversing has become an important method widely us ed in surveys. A traverse is represented
by a series of straight lines called traverse legs, which are used to connect a set of selected
points on the earth, called traverse stations. The surveyor makes length and direction
measurements and uses them to compute the relative positions of the traverse stations in a
coordinate system. Meantime, the surveyor must make the adjustment work of the indirect
observations that compose the traverse and for this operation the least square method of
indirect observations is used. During the entire pr ocess of the network adjustment, computing
the approximate points coordinates is one of the im portant steps [2].
For the present study, the threeAdimensional geodet ic network is represented by a
closedAloop traverse, which starts at a point and e nds at the same point station. In order to
minimize the effect of the accidental errors made i n measurements, they are all adjusted by
several computations. For the present situation, th e type of closedAloop traverse is preferred,
because it provides a check over the fieldwork, com putations and starting data. Meantime, it
provides a basis for comparison, in order to determ ine the overall accuracy of the entire work.
Since the main purpose of a traverse is to locate p oints relative to each other on a
common grid, there are needed certain elements of s tarting data. For this study, the initial data
is represented by the coordinates of the starting/e nding point station.
The present work develops an adjustment application of the geodetic measurements,
designed with the help of the Matlab programming la nguage. It is presented the mathematical
model of least squares method and all the phases in achieving the results of the proposed
algorithm. For a data set consisting in three diffe rent types of observations, like: length
observations, direction observations and angles, th ere were determined the most probable
values for the points coordinates with the accuracy of their determination and the most
probable values for the measurements, of course wit h their associated accuracy [3].
The current work presents the steps taken to create a program that automatically
processes the measurements from a threeAdimensional geodetic network, in which points are
defined by all the three coordinates (3D Space).

2. Presentation of the mathematical algorithm

2.1. Starting data

For the present work the geodetic network is repre sented by a closedAloop traverse,
located approximately in the center of Iasi city (R omania), surrounding the “Dosoftei House”
Museum of Old Moldavian Literature (Fig. 1).

Fig. 1. The closed A loop traverse scheme

A. M. Loghin, R.S Ajin, V. E. Oniga
The threeAdimensional geodetic networks adjustement automation using Matlab

17
The traverse has two points with known coordinates (500 and 501) and also three
unknown points (502, 503, 504) that must be compute d. The known points coordinates
resulted from a previous stage, using GNSS real tim e kinematic measurements, made with a
South S82 – V receptor. During the traverse, the ho rizontal angles, vertical angles, and
distances are measured, using a total station, with an accuracy of 7'' for angles respectively of
2 mm + 2 ppm for distances.

2.2. Data processing

In order to create the adjustment application, the Matlab programming language was
used. Matlab is a highAperformance technical comput ing environment, developed by
MathWorks, which combines both math and graphic fun ctions with a powerful highAlevel
language.
Some of the most important features of this languag e are: efficient matrix and vector
computations, easy creation of scientific and engin eering graphics, plotting of functions and
data, algorithms implementation, application develo pment including graphical user interface
building, objectApriented programming and interfaci ng with programs written in other
languages, including C, C++, Java, Fortran and Pyth on.
Matlab can be used for a range of applications, inc luding signal processing and
communications, image and video processing, control systems, test and measurement,
computational finance and computational biology. Be cause of its important feature of matrix
and vector computation and of algorithm manipulatio n, it is primarily used in producing
solutions to complex systems of equations, in model ing, simulation and prototyping and also
in data analysis, exploration and visualization. Mo re than a million engineering and scientists
in industry and academia use Matlab, the language o f technical computing [4].

2.3. Mathematical support

In order to compute the most probable values for po ints coordinates with the accuracy
of their determination and also the corrections app lied to the observed angles and measured
distances, the least square method of indirect obse rvations was used. This is the most rigorous
adjustment procedure, based on the mathematical the ory of probability that the sum of the
squares of the errors times their respective weight s are minimized. When vi is assumed as a
small correction value to be added to each observed value and pi as a weight of observed
values, each value of vi is determined in a manner to satisfy the following formula [5]:

2
i i
ip v ∑ = minimum (1)
In this program, one observation equation is develo ped for one set of observations and
after theoretically evaluating the weight of observ ations, each equation is solved by the least
squares method in order to satisfy equation (3) [6] .
An important aspect of the least square method is t hat it permits to all kind of
observations, regardless of their number or type, t o be included into the adjustment process
and simultaneously used in the computations. Thus, the adjustment combines distances,
horizontal angles, azimuths, zenith or vertical ang les, coordinates and also GNSS
observations.
The entire adjustment process is composed by severa l steps, such as: computing the
provisional values of the coordinates, preparing an d solving the functional model of
corrections, computing the adjusted elements and fi nally the accuracy analysis.

„1 Decembrie 1918” University of Alba Iulia RevCA D 19/2015

18
The functionalAstochastic model for the indirect we ighted measurement method,
containing a number of r equations with n unknown parameters is represented by the
following equation [5]:
BrAn XnA1 + L rA1 = V rA1, with P rAr weighting matrix (2)

where: B rAn is the coefficient matrix of the correction equati ons,
X nA1 is the vector matrix of the unknown parameters (th e corrections of the unknown
point coordinates in a local geodetic system),
L rA1 is the vectorAmatrix of the free elements,
V rA1 is the correction matrix.

All the correction coefficients and also the free t erms were calculated for the input
data set composed by: distances, horizontal and ver tical angle measurements.
The detailed algorithm of the threeAdimensional geo detic network adjustment with its
mathematical implementation can be entirely found i n [7].

3. Results and discussion

3.1. Logical support

Before starting the implementation process of the a djustment in the Matlab
programming language, the mathematical algorithm wa s carefully analyzed and also it was
realized the flow chart of the program. A rough flo w of the program is presented in Fig. 2.

Fig. 2. Rough flow of the program

As input data, there was used two text files: one w ith the known point coordinates and
the second one containing the measurements acquired with the help of a total station.

A. M. Loghin, R.S Ajin, V. E. Oniga
The threeAdimensional geodetic networks adjustement automation using Matlab

19

(a) (b)

Fig. 3. Model of the input files, (a) point coordin ates text file ( Points.txt ), (b) Measurements
text file ( Traverse.txt )

The present program retains all measured data from the text file in a cellAarray D,
having on its columns the observations in the follo wing order: the name of the points,
distances, vertical angles, horizontal directions, reflector elevation, code. In the total stations
output text files, all observations are organized b y point stations, being delimited by the
“99999 0 0 0 0 0 0” sequence. Therefore, in order t o store the data from each point station, the
following code lines sequence was implemented:

fid = fopen( 'Traverse.txt' ,'rt' );
C = cell2mat(textscan(fid, '%f%f%f%f%f%f' ));
fclose(fid);
idx = [1,diff(C(:,1)'~=9999)];
idy = find(idx<0)Afind(idx>0);
idy(2,:) = 1;
D = mat2cell(C,idy(:),size(C,2));
D(2:2:end) = [];

Fig. 4. Extracting the observations from each point station

Next, there were defined the start, intermediate an d the end point station points
respectively, in order to determine the traverse di rection. The dialog box allows the user to
input the points’s names from the keyboard (Fig. 5) .

Fig. 5. Input dialog box

In the next steps, the application computes the pro visional coordinates for all the
unknown points of the traverse. Then, following the mathematical algorithm, the program
realizes the adjustment process of the threeAdimens ional geodetic network. For each stage,
there were created new script files, that were impl emented as functions in the main program.
Therefore, there were created several functions suc h as: the calculation of the horizontal
angles function, the measurement errors calculation function, conversion function between
geodetic and rectangular coordinates, azimuth angle function, coefficient calculation
functions, accuracy analysis functions.

„1 Decembrie 1918” University of Alba Iulia RevCA D 19/2015

20
3.2. Conversion function between geodetic and recta ngular coordinates

Based on the geodetic coordinates (B, L) and the co nstant coefficients a 00 , a 10 , a 20 , a 30 ,
…, b 01 , b 11 , b 21 , there were computed the plane rectangular coordin ates in the National
Projection System, called „Stereographical on uniqu e secant plan” and the normal altitudes in
the „Black Sea 1975”. The conversion algorithm is d escribed in [8].
Firstly, using the particular function “dms_degrees ” there were computed the geodetic
coordinates from degrees, minutes, seconds into deg rees with their decimal part, the results
being stored in the ‘Points_BLH.txt' text file.

function [ degreesB,degreesL ] = dms_degrees( dmsB,dmsL )
disp( ‘Degrees, minutes seconds to degrees’ )
%input data
Coord =dlmread( 'Points.txt' );
dmsB=Coord (:,2:4);
dmsL=Coord (:,5:7);
He=Coord (:,8);
%end of input data
format long g
degreesB = dms2degrees(dmsB);
degreesL = dms2degrees(dmsL);
coord_BLH(:,1)=Coord_cun(:,1);
coord_BLH(:,2)=degreesB;
coord_BLH(:,3)=degreesL;
coord_BLH(:,4)=Coord_cun(:,8);
dlmwrite( ‘Points_BLH.txt' ,coord_BLH, 'delimiter' ,'\t' ,'newline' ,'pc' ,'precision, 9);
end

Fig. 6. The conversion function from degrees, minut es, seconds into degrees

The Matlab function uses the constant coefficients, stored in two different arrays
par_a and par_b . Further, there are defined the geographical coord inates of the central point
of Romania (φ 0, λ 0) and there were created the B and L vectors that contain the points
latitudes and longitudes respectively. Also, there were computed the differences between
points coordinates and the central point coordinate s, that were stored in the f and l vectors.
In the next stage, there were computed the addition al arrays Fa , Fb , La , Lb and finally,
the Xtg and Ytg vectors that contain the rectangular coordinates o n the tangent plane were
calculated, using the following equations [8]:

X TG = ( a 00 + a 10 f + a 20 f2 + a 30 f3 + a 40 f4 + a 50 f5 +a 60 f6 ) 1.000 +
( a 02 + a 12 f + a 22 f2 + a 32 f3 + a 42 f4 )l 2 +
( a 40 + a 14 f + a 24 f2 ) l 4 +
( a 06 ) l 6 +……… (3)
Y TG = ( b 01 + b 11 f + b 21 f2 + b 31 f3 + b 41 f4 + b 51 f5 ) 1,000 +
( b 03 + b 13 f + b 23 f2 + b 33 f3 ) l 3 +
( b 05 + b 15 f ) l 5 + ……….

These coordinates are passed into the secant plane of the projection, by multiplication
with the coefficient of reduction to scale (C = 0.9 99750) and adding the translation of the
coordinates system axes (X = X tg C+X o; Y = Y tg C+Y o).

A. M. Loghin, R.S Ajin, V. E. Oniga
The threeAdimensional geodetic networks adjustement automation using Matlab

21

function [ coord_stereo ] = conversionBLH_XYZ (Coord)
par_a = dlmread( 'parametri_a.txt' );
par_b = dlmread( 'parametri_b.txt' );
pct = Coord; j = size(pct,1);
B0=46; L0=25;
B=pct(:,2); L=pct(:,3);
for i=1:j
f(i,1)=3600*(B(i,1)AB0)*10^A4; l(i,1)=3600*(L( i,1)AL0)*10^A4;
end
for i=1:j
F(:,i)=[f(i,1)^0;f(i,1)^1;f(i,1)^2;f(i,1)^3;f(i, 1)^4;f(i,1)^5;f(i,1)^6;f(i,1)^7];
Fb(:,i)=[f(i,1)^0;f(i,1)^1;f(i,1)^2;f(i,1)^3;f(i ,1)^4;f(i,1)^5;f(i,1)^6];
La(:,i)=[l(i,1)^0;l(i,1)^2;l(i,1)^4;l(i,1)^6;0;0 ;0;0];
Lb(:,i)=[l(i,1)^1,l(i,1)^3,l(i,1)^5,l(i,1)^7,0,0 ,0];
end
Pa=par_a*F; Xtg=dot(Pa,La);
Xsec=Ntg*0.999750; X=500000+Xsec;
Pb=par_b*Fb; Ytg=dot(Pb,Lb);
Ysec=Ytg*0.999750; Y=500000+Ysec;
STEREO (:,1)=pct(:,1); STEREO (:,2)=X; STEREO (:,3) =Y;STEREO (:,4)=pct(:,4);
dlmwrite( 'Points_stereo.txt' ,STEREO, 'delimiter' ,'\t' ,'newline' ,'pc' ,'precision' ,9);
Points_stereo=[STEREO (:,1) STEREO (:,2) STEREO (:, 3) STEREO (:,4)];
end

Fig. 7. Conversion function between geodetic and re ctangular coordinates

3.3. Population of the weighting matrix procedure

P=zeros(n,n);
for i=1:length(d)
for j=1:length(d)
if i==j
P(i,j)=weight_s; k=i; end ; end
end
for i=(k+1):(k+length(weight _beta))
for j=(k+1):(k+length(weight_beta))
if i==j
P(i,j)= weight _beta(iAn); k=i; end ; end
end
for i=(k+1):n
for j=(k+1):n
if i==j
P(i,j)= weight _z; end ; end
end
Fig. 8. Population of the weighting matrix procedur e

The weighting matrix is a diagonal matrix that cont ains the weight values for each
measured data. The procedure contains three for ins tructions, in order to select the weight
values for distances, horizontal and vertical angle s respectively.

„1 Decembrie 1918” University of Alba Iulia RevCA D 19/2015

22

4. Conclusions

In recent years automation, computer science and a pplied informatics domains have
known a great evolution, being used in almost all f ields of activity, such as: engineering,
archaeology, surveying, medical and chemical indust ries, design projects.
The aim of this study was to obtain an automated m ethod for processing the
measurements from a threeAdimensional geodetic netw ork, in the form of a closedAloop
traverse, using Matlab programming language.
In order to prepare the entire program, both knowle dge of geodesy (threeAdimensional
networks, processing of the measurements from a thr eeAdimensional geodetic network) and
programming (Matlab programming language) are requi red. Firstly, the mathematical
algorithm of adjustment was carefully analyzed and there was created the flowchart of the
program.
After processing, the program offers the solutions: the most probable values for the
points coordinates with their accuracy determinatio n and the most probable values for the
measurements also with their associated accuracy.
This program has more advantages, such as: it has a friendly interface being easy to
understand, utilizes simple notations which allows focusing on mathematical algorithm, has a
great computing speed, being composed by more parti cular functions and procedures.
In the execution process of the program, there wer e achieved important properties,
like: correctness, compatibility, portability, effi ciency, clarity.
Programming represents an essential domain of info rmatics. It helps us to understand
the hipperAconnecting world we are living in and to represent what happens in the backside of
the displays, in a coded language.

5. References

1. Ghilani D.C. V Adjustement computations. Spatial da ta analysis. Fifth edition, John Wiley
& Sons, Inc., Hoboken, ISBN 978V0V470V46491V5, New Jersey, S.U.A, 2010
2. Shi LiVmei, Chen YiVjin, Gong Hui, Wang Ning – Algo rithm design and realisation of the
approximate coordinate of points in the computation of traverse network indirection
adjustment, University of Mining and Technology, En gineering of surveying and mapping,
Beijing 100083, China, June 2006
3. A. Călin, A. Gînța, P. D. Dumitru – Programming the free adjustment of the precise
leveling measurements, Journal of Geodesy and Cadas tre RevCAD, no. 16, pag. 26V33,
AlbaVIulia, Romania, ISSN 1583V2279, 2014
4. Matlab user guide. Source: https://www.mathworks.com/help/pdf_doc/matlab/getst art.pdf
5. Nistor Gh., Teoria prelucrării măsurătorilor geodez ice, Universitatea Tehnică “Gh.
Asachi”, Iași, 1995
6. Takehisa Harada – Universal program for adjustment of any geodetic network (PAGVU),
Journal of the Geodetic Society of Japan, Vol. 26, No. 3, pp.147V171, 1980
7. C. Chirilă – “Geodezie matematică – Îndrumător de l ucrări practice și proiect”,
Publishing house Tehnopress, Iași, Romania, 2014
8. Moca Valeriu, Chirilă Constantin V Cartografie mate matică. Întocmire și redactare hărți ,
Lecture Notes, Technical University „Gheorghe Asach i” of Iasi, Faculty of Hydrotechnical
Engineering, Geodesy and Environmental Engineering, 2002