3D Geodetic Least Squares Adjustment Programs

[Eugene Ballance]

 For several years, I have also been working with a 3D program based on C. Ghilani's book. However, mine uses Octave(Matlab) scripts. Yours is the first open source code that I could  compare my scripts with. Really happy when I found you!.

My observations are initially entered into shapefiles in Christine-GIS (closest equivalent to ArcView3.x that works with Windows 10) which point to hand written copy images of  original US Coast Survey Triangulation Documents. 

Octave scripts read the ".dbf" portion data, from the shapefiles, for the subset of the records that are to be used in a given adjustment. Octave is good with matrix operations, and I am interested to find how it compares with Eigen.

I have attached my most current working version. I have now included a process to convert input data to a SALSA project, which can be opened in the SALSA gui.

I am still not clear on how the "n","e", and "u" components can be fixed individually. The seems to be no example of this in the user's manual. As you may know, in the 3d case Mr. Ghilani exclusively used small SDs to do this.

An example of this would be very helpful.

Keep up the good work!

~Gene Ballance

[Eugene Ballance]

I was able to fix the "u" components. I will soon post a Google Drive link  for the SALSA file results for the Pamlico Sound triangulation & oyster studies data.

In my Octave program, the sequence converges, but I used small SDs, instead of fixing n/e/u components.

Seems like SALSA may have given up too quickly? I guess I can change the number of iterations.

I also did some study of the NGS program "ADJUST.exe". As of my checking last year, there was no way to use it unless one used at least one pseudo GPS observation, in addition to terrestrial. I had no easy way to debug it at the time, so I started on writing the Octave version I am using now.

On this problem, the Octave version seems to run much faster than SALSA or Ghilani's "Adjust"..if you think this matters. 

~Gene Ballance

[Eugene Ballance]

Link to "fairfield Pamlico Sound" files...

https://drive.google.com/drive/folders/1jderEgxyK8O5khk0xMoinheLVdA3emSE?usp=sharing

[Eugene Ballance]

The data in this project was gathered instruments with widely different accuracy..10" theodolites( 5 SOA) with many repetitions and averaging, marine sextants( 80 SOA) with little redundancy, and ordinary theodolites whose accuracy may not have been much more than the sextants. Thus if I set in the "project.cfg" file a convergence interval of 1.0e+0, then I get "convergence". 

Also, the sextant fixes were mostly resections with just two angles measured and often close to the "danger circle", so one might wonder why I would try to include them in the same adjustment. The reason is that I have found that many of the shore signals were also occupied with the sextants, so that including all the fixes tends to "harden" the triangulation network, if more than three signals were observed when doing a resection fix.

Latest adjustment at...

https://drive.google.com/drive/folders/15aQumHqqD3f0SsFqwpU3QN3k03LmLdFd?usp=sharing

 

[Eugene Ballance]

Another observation that the oyster surveys made frequently was a transit, or zero angle observation of two signals. I tried to add this in my 2D Octave adjustment script some years ago, but it seemed to be very unstable. I am now trying it in my Octave 3D and SALSA and there seems to be no problem, if I add  a pseudo ZANG. 

---Has anyone worked with these "transits" in least squares adjustment programs?

My strategy is to run the transit with a 90 deg observation ZANG on the first run, and then adjust it as to the residual on subsequent  runs to minimize the residual. I hope this will add enough stability to the network.

[Eugene Ballance]

Ha! Seems there is still a problem with the zero angles.. on the larger networks,anyway.

There is not a problem on the simple networks like The Ghilani 3D example.

But on the large Pamlico Sound network even using pseudo SDs of 0.001 on the HANGs doesn't help.

As I now recall, I ended up just calling these "external" (to the adjustment) and recalculating them after an adjustment.

They are usually just one observation, anyway, not many that need to be adjusted to a most probable value.

[Eugene Ballance]

If careful to put initial value on the correct side of the transit line, divergence is only slowed a bit.

I have often entertained the possibility of adding these as constraints on the XYZ coordinates, so that for example, X & Z would have their own slopes with respect to Y, then take the differentials. Going to NEU would be via the usual rotation matrix. 

This would basically add another type of observation.

[Corwin Salsa]

Hi Eugene,

Thanks for reaching out! Looks like you’ve been working today through many of your original questions for how to employ SALSA for your project. I will admit that I do not follow all of the details, but below are a few ideas.

Regarding how SALSA fixes individual north, east, or up components, it looks like you’ve figured out how to do that in the UI, but section 4.2.4 from the SALSA manual describes how these constraints are performed: “These constraints are implemented in a strictly rigorous way using matrix decomposition and projection. This has the effect of limiting the vector space of the least squares solution to the subspace in which the constraints are exactly satisfied. This is unlike other techniques for applying constraints that use “overweighting” or pseudodata with zero covariance to force the least squares solution to hold components fixed [1]; these methods are inexact and run the risk of destabilizing the problem and even making it singular.” We anticipate providing more mathematical detail in a future SALSA release. 

Regarding the convergence criteria, it looks you’ve determined how to change both the convergence threshold and the maximum number of iterations, but let us know if you continue to have problems.

Regarding your statement: “The data in this project was gathered instruments with widely different accuracy..10" theodolites( 5 SOA) with many repetitions and averaging, marine sextants( 80 SOA) with little redundancy, and ordinary theodolites whose accuracy may not have been much more than the sextants. Thus if I set in the "project.cfg" file a convergence interval of 1.0e+0, then I get "convergence".”

I have not had a chance to look through the SALSA project you’ve posted, but remember that SALSA provides the ability to specify different measurement uncertainties for different measurements and measurement types. E.g. you can create a specific UNCR record for each measurement type and then tie those to the corresponding measurements. That may allow you to converge on your solution more correctly and without the need to modify the convergence criteria. Of course if you’ve already done this, then my apologies for the redundancy. 

I do not know if these “transits” you’ve described have been studied in detail, but in the real world I would be surprised if you ever got a true zero angle observation - just very small. Unless of course I do not understand the measurement you’re describing. I would also be cautious about adding pseudo-measurements to make a LS process work - that can easily lead you astray. Usually when pseudo-measurements are required as you describe, there are fundamental challenges in the underlying geometry that must be addressed (e.g. via constraints, etc.) Sounds like from your last post that you were thinking along the same lines. And regarding your comment “This would basically add another type of observation.”, I think it’s more like you’re adding known information (in the form of the constraint), and thus improving the solution - but it’s not really an observation in the true sense of the word. But I suspect we’re on the same page there as well, from your context.

Good luck, and let us know if you have further questions - though we may have to do a bit of translation of the technical details to assess more detailed problems. This kind of dialog is exactly why we set up this group.

Thanks!

Corwin

[Eugene Ballance]

Thanks so much for your excellent response!

Of course, no measurement is exact. Maybe what I am envisioning  is a constraint on parameters. On the other hand, aren't all the other observations written in terms of parameters? Just saying. The paper I got this idea from is attached. What I should have said is:

  "There will be another 9D function involving  parameters to take partial derivatives of, in the adjustment code".

Problem is, when I write this in the code, it is indistinguishable from a "zero" observation. That is, the procedure Iooks exactly like the one for the other observations ..DIST, AZIM, HANG  in that it is "thrown into the mix". It will contribute to, but not force, the solution. I will give it a an initial SD (which will have to be in m^2). The question then becomes "how did you measure that". My answer is "I don't know, it might be the area swept out by a line centered on the observation position (i) and connected to the (b) and (f) positions, given their probable errors". It would also include an instrument centering error, so it could be modeled as in the HANG case, but with areas (m^2) calculated instead of SOA.

My answer to Mr. Harvey's question on page 9 (of the attached) is.."add the equations together". This is what I am now trying to implement in my Octave version.

I started out with FORTRAN in the early 1970's. The computer languages used since then have been HP48 RPL, Avenue (in Arcview 3.x), the scripting language of Christine-GIS (something like c++, have to declare variables), and Octave(close to Matlab). With some work, I could probably use Octave exclusively. Of all of them, this is my preference, because it is based on matrix (vector) operations. I will probably get into c++ more now, since you have so graciously posted your code (btw..is all the code posed? ).

Thanks for your guidance!

~Gene

[Eugene Ballance]

Still working on a way to calculate an apriori SD for my slope "observations". The function is ...

f(xb,yb,zb,xi,yi,zi,xf,yf,zf) = (((yf-yi)*(xi-xb))-((yi-yb)*(xf-xi)))+(((yf-yi)*(zi-zb))-((yi-yb)*(zf-zi))) 

..which = 0 when you are on the 3D line-of-sight of two stations (or you would be from one of the other two).

Partials are pretty simple. I have a function called "getvfact.m" that will give fractions of distance from two pairs of points&3-vectors to point of closest approach of the lines, so I can use this with direct3d.m and inverse 3d.m to get an estimated position.

Plan is then to use the usual sextant SD and the est. position to get a max XYZ offset. Then plug the est. pos into f & compare the value obtained with that of just adding the offset to xi,yi,zi and calculating the f-value.

~Gene

[Eugene Ballance]

Now ready for the slope observation test, but it will have to wait until I have time to debug it, if it fails.

[Eugene Ballance]

-failed big-time. Seems very unstable. 

Now revisiting my previous method of handling the zero-angle observations, which was to add angle observations of 59 SOA from bs-i-fs AND fs-i-bs. This did converge in the large networks.

I also want to add the timed observation points on transects between two fixes. The equations needed for this are a pair of the near-zero-angle-observations, the equality of two zenith angle observations (zi = zj, so zi-zj = 0), and that a ratio of 2 distances is a given constant (di/dj = k, so di-k*dj = 0).  I am looking into including these as constraints as P.A. CROSS did in his paper, using the method of Lagrange multipliers. The result being that I would calculate probable error ellipses for the positions of timed observation points.

~Gene

[Eugene Ballance]

So glad to have SALSA to compare my Octave results. Found one major blunder (6 DOA!) in one the 1886 angle records. This helped a lot!  I am using a large std on blunders, instead of eliminating them. Attached is my latest effort.

~Gene