A Multifidelity Approach to Robust Orbit Determination

Alberto Foss\`a; Roberto Armellin; Emmanuel Delande; Matteo Losacco; and Francesco Sanfedino

arXiv:2302.14655·math.NA·November 7, 2023

A Multifidelity Approach to Robust Orbit Determination

Alberto Foss\`a, Roberto Armellin, Emmanuel Delande, Matteo Losacco, and Francesco Sanfedino

PDF

Open Access

TL;DR

This paper introduces a multifidelity algorithm that enhances orbit determination robustness by refining initial estimates, detecting outliers, and efficiently propagating uncertainties using differential algebra and domain splitting techniques.

Contribution

The paper presents a novel multifidelity approach combining differential algebra and automatic domain splitting to improve orbit determination robustness and outlier detection.

Findings

01

Effective outlier detection in orbit data

02

Improved initial orbit estimates with reduced uncertainty

03

Validated on synthetic and real TAROT network data

Abstract

This paper presents an algorithm for the preprocessing of observation data aimed at improving the robustness of orbit determination tools. Two objectives are fulfilled: obtain a refined solution to the initial orbit determination problem and detect possible outliers in the processed measurements. The uncertainty on the initial estimate is propagated forward in time and progressively reduced by exploiting sensor data available in said propagation window. Differential algebra techniques and a novel automatic domain splitting algorithm for second-order Taylor expansions are used to efficiently propagate uncertainties over time. A multifidelity approach is employed to minimize the computational effort while retaining the accuracy of the propagated estimate. At each observation epoch, a polynomial map is obtained by projecting the propagated states onto the observable space. Domains that do…

Tables22

Table 1. Table 1 : Epoch, duration and number of measurements for the considered passages

Sensor	Epoch (UTC)	Duration	#measurements
La Reunion	2019-02-25T18:49:01.148	05:13:10.189	8
La Reunion	2019-02-27T23:27:20.511	00:00:48.972	3
La Reunion	2019-03-01T01:35:55.561	00:00:48.019	3
La Reunion	2019-03-01T23:10:59.277	00:00:24.011	2
Calern	2019-03-01T23:58:27.006	00:00:24.926	2

Table 2. Table 2 : Root mean square error for a five days propagation

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
LF	$181.501 633 403 219$	$121.951 277 883 348$	$6.504 145 374 968 15$	$0.031 901 313 098 891 3$	$0.010 522 041 666 980 1$	$0.000 184 099 443 676 873$
MF	$1.810 795 649 502 28$	$1.321 405 755 533 21$	$0.085 518 923 696 812 7$	$0.000 301 435 053 817 822$	$0.000 101 002 771 932 697$	$2.811 987 724 767 23 × 10^{- 06}$
HF	$0.140 434 694 853 268$	$0.101 660 903 665 224$	$0.006 797 315 686 115 94$	$2.290 873 391 004 1 × 10^{- 05}$	$7.563 036 510 961 41 × 10^{- 06}$	$4.013 605 887 364 88 × 10^{- 07}$

Table 3. Table 3 : Computational time for a five days propagation

	Low-fidelity	Multifidelity	High-fidelity	Single sample	Monte Carlo
$t$ , $s$	$0.237$	$3.392$	$144.934$	$3.040 987 5$	$30 409.875$
$t / t_{M F}$	$0.069 870 283 018 867 9$	$1$	$42.728 183 962 264 2$	$0.896 517 541 273 585$	$8965.175 412 735 85$

Table 4. Table 4 : Initial conditions for target and outlier orbits

	$a$ , $km$	$e$ , -	$i$ , $°$	$ω$ , $°$	$Ω$ , $°$	$M_{0}$ , $°$
Target	$22 953.852 669 768 778$	$0.707 854 612 716$	$3.387 521 317 683$	$172.980 213 527 756$	$- 168.891 499 315 499$	$60.742 995 057 860$
Outlier	$22 953.852 669 768 778$	$0.687 854 612 716$	$3.387 521 317 683$	$172.980 213 527 756$	$- 168.891 499 315 499$	$60.742 995 057 860$

Table 5. Table 5 : Test case scenarios with simulated measurements

Scenario #	Pruning	Outliers
A	NO	NO
B	YES	NO
C	NO	3^rd passage
D	YES	3^rd passage

Table 6. Table 6 : Number of domains as function of propagation time for scenario B

	Number of domains
Time since $t_{0}$ , $h$	#1 Propagation	#2 Projection	#3 Pruning	#4 Merging
$0$	1	1	1	1
$0.006 673 888 888 888 89$	1	1	1	1
$0.013 325 833 333 333 3$	1	1	1	1
$4.393 378 888 888 89$	1	1	1	1
$4.400 048 333 333 33$	1	1	1	1
$5.206 175 277 777 78$	1	1	1	1
$5.212 832 222 222 22$	1	1	1	1
$5.219 496 944 444 44$	1	1	1	1
$52.638 711 944 444 4$	1	9	7	5
$52.645 372 777 777 8$	5	7	7	5
$52.652 315 277 777 8$	5	7	7	5
$78.781 781 388 888 9$	5	63	38	18
$78.788 455 277 777 8$	18	38	38	18
$78.795 12$	18	38	38	18
$100.366 146 944 444$	18	38	38	18
$100.372 816 666 667$	18	38	38	18
$101.157 182 777 778$	18	38	35	19
$101.164 106 666 667$	19	35	35	19

Table 7. Table 7 : Number of domains as function of propagation time for scenario D

	Number of domains
Time since $t_{0}$ , $h$	#1 Propagation	#2 Projection	#3 Pruning	#4 Merging
$0$	1	1	1	1
$0.006 673 888 888 888 89$	1	1	1	1
$0.013 325 833 333 333 3$	1	1	1	1
$4.393 378 888 888 89$	1	1	1	1
$4.400 048 333 333 33$	1	1	1	1
$5.206 175 277 777 78$	1	1	1	1
$5.212 832 222 222 22$	1	1	1	1
$5.219 496 944 444 44$	1	1	1	1
$52.638 711 944 444 4$	1	9	7	5
$52.645 372 777 777 8$	5	7	7	5
$52.652 315 277 777 8$	5	7	7	5
$78.781 781 388 888 9$	5	63	63	5
$78.788 455 277 777 8$	5	63	63	5
$78.795 12$	5	63	63	5
$100.366 146 944 444$	5	63	38	18
$100.372 816 666 667$	18	38	38	18
$101.157 182 777 778$	18	38	35	19
$101.164 106 666 667$	19	35	35	19

Table 8. Table 8 : Success rate of outlier detection with simulated measurements

Case #	Dynamics	True positives	False negatives	False positives	True negatives
A,B	LF	13	5	0	0
A,B	MF	18	0	0	0
C,D	LF	10	5	0	3
C,D	MF	15	0	0	3

Table 9. Table 9 : True state and initial guess for scenarios A to D

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
True	$- 21 551.184 664 630 193$	$14 404.866 452 074 804$	$- 1082.462 558 770 526$	$- 3.580 403 901 491$	$- 0.736 464 589 895$	$0.001 943 794 765$
Guess	$- 21 553.629 554 466 817$	$14 406.082 502 818 883$	$- 1082.362 628 418 208$	$- 3.579 945 419 480$	$- 0.736 567 402 660$	$0.001 921 790 734$

Table 10. Table 10 : LS and LSAR estimation errors and ± 3 σ plus-or-minus 3 𝜎 \pm 3\sigma bounds for scenario A and B

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
LS Error	$0.001 432 882 666 337 62$	$0.166 254 060 672 003$	$0.064 043 989 450 272 4$	$1.436 626 695 255 73 × 10^{- 05}$	$6.952 457 072 228 16 × 10^{- 06}$	$8.227 634 197 685 02 × 10^{- 07}$
$\pm 3 σ$ bounds	$0.911 943 829 770 987$	$0.456 451 793 548 531$	$0.205 293 174 909 563$	$0.000 171 046 554 685 692$	$4.022 577 098 770 14 × 10^{- 05}$	$1.534 972 270 643 65 × 10^{- 05}$
LSAR Error	$0.180 380 981 372 442$	$0.082 141 679 550 848$	$0.045 135 145 490 439$	$1.904 610 809 208 85 × 10^{- 05}$	$1.059 628 500 081 36 × 10^{- 05}$	$2.082 979 514 912 71 × 10^{- 07}$
$\pm 3 σ$ bounds	$1.130 409 923 998 03$	$0.565 739 521 827 316$	$0.254 563 548 634 605$	$0.000 212 009 774 878 624$	$4.987 191 556 960 88 × 10^{- 05}$	$1.903 318 934 209 41 × 10^{- 05}$

Table 11. Table 11 : LS and LSAR estimation errors and ± 3 σ plus-or-minus 3 𝜎 \pm 3\sigma bounds for scenario C

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
LS Error	$579.080 210 298 213$	$377.552 242 082 303$	$34.482 667 119 512 7$	$0.118 092 221 182 537$	$0.024 624 382 030 599 6$	$0.002 568 627 842 706 61$
$\pm 3 σ$ bounds	$0.861 941 148 124 937$	$0.422 408 271 524 913$	$0.210 029 303 800 098$	$0.000 157 272 736 644 034$	$3.680 251 141 795 16 × 10^{- 05}$	$1.472 566 609 882 91 × 10^{- 05}$
LSAR Error	$0.772 982 128 250 064$	$0.257 594 372 595 664$	$0.014 885 077 288 853 4$	$0.000 138 442 689 863 824$	$4.156 257 135 901 47 × 10^{- 05}$	$3.410 044 251 506 9 × 10^{- 07}$
$\pm 3 σ$ bounds	$1.130 151 010 836 91$	$0.565 569 011 378 571$	$0.254 564 917 720 344$	$0.000 211 949 556 622 479$	$4.986 384 705 859 85 × 10^{- 05}$	$1.903 236 267 704 × 10^{- 05}$

Table 12. Table 12 : LS and LSAR estimation errors and ± 3 σ plus-or-minus 3 𝜎 \pm 3\sigma bounds for scenario D

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
LS Error	$0.188 526 166 259 416$	$0.043 656 216 609 869$	$0.044 030 023 872 879 9$	$2.445 104 726 108 99 × 10^{- 05}$	$1.512 146 090 805 39 × 10^{- 05}$	$2.041 683 136 327 85 × 10^{- 06}$
$\pm 3 σ$ bounds	$1.096 327 622 030 42$	$0.598 195 189 390 871$	$0.226 457 257 104 369$	$0.000 211 436 830 082 596$	$4.794 845 777 630 55 × 10^{- 05}$	$1.678 986 243 955 95 × 10^{- 05}$
LSAR Error	$0.333 753 328 664 473$	$0.005 784 302 300 751 49$	$0.055 575 712 109 970 7$	$4.972 276 519 616 83 × 10^{- 05}$	$1.859 398 898 995 53 × 10^{- 05}$	$4.301 013 176 568 8 × 10^{- 07}$
$\pm 3 σ$ bounds	$1.359 341 110 884 77$	$0.741 705 821 914 885$	$0.280 806 887 552 135$	$0.000 262 159 299 521 155$	$5.945 448 104 589 99 × 10^{- 05}$	$2.081 914 492 261 71 × 10^{- 05}$

Table 13. Table 13 : Runtime and number of iterations with simulated measurements

Case #	Pruning	LS	LSAR	LS	LSAR
	$t$ , $s$			# iterations
A	-	$331.642$	$125.777$	3	2
B	$70.297$	$331.257$	$125.944$	3	2
C	-	$642.552$	$324.101$	6	5
D	$40.718$	$133.349$	$87.541$	2	2

Table 14. Table 14 : Test case scenarios with real measurements from the TAROT network

Case #	Pruning	Outliers
E	NO	NO
F	YES	NO
G	NO	3^rd passage
H	YES	3^rd passage

Table 15. Table 15 : Number of domains as function of propagation time for scenario F

	Number of domains
Time since $t_{0}$ , $h$	#1 Propagation	#2 Projection	#3 Pruning	#4 Merging
$0$	1	1	1	1
$0.006 673 888 888 888 89$	1	1	1	1
$0.013 325 833 333 333 3$	1	1	1	1
$4.393 378 888 888 89$	1	1	1	1
$4.400 048 333 333 33$	1	1	1	1
$5.206 175 277 777 78$	1	1	1	1
$5.212 832 222 222 22$	1	1	1	1
$5.219 496 944 444 44$	1	1	1	1
$52.638 711 944 444 4$	1	3	3	1
$52.645 372 777 777 8$	1	3	3	1
$52.652 315 277 777 8$	1	3	3	1
$78.781 781 388 888 9$	1	81	45	15
$78.788 455 277 777 8$	15	45	45	15
$78.795 12$	15	45	45	15
$100.366 146 944 444$	15	135	116	24
$100.372 816 666 667$	24	116	116	24
$101.157 182 777 778$	24	116	116	24
$101.164 106 666 667$	24	116	116	24

Table 16. Table 16 : Number of domains as function of propagation time for scenario H

	Number of domains
Time since $t_{0}$ , $h$	#1 Propagation	#2 Projection	#3 Pruning	#4 Merging
$0$	1	1	1	1
$0.006 673 888 888 888 89$	1	1	1	1
$0.013 325 833 333 333 3$	1	1	1	1
$4.393 378 888 888 89$	1	1	1	1
$4.400 048 333 333 33$	1	1	1	1
$5.206 175 277 777 78$	1	1	1	1
$5.212 832 222 222 22$	1	1	1	1
$5.219 496 944 444 44$	1	1	1	1
$52.638 711 944 444 4$	1	3	3	1
$52.645 372 777 777 8$	1	3	3	1
$52.652 315 277 777 8$	1	3	3	1
$78.781 781 388 888 9$	1	81	81	1
$78.788 455 277 777 8$	1	81	81	1
$78.795 12$	1	81	81	1
$100.366 146 944 444$	1	243	116	24
$100.372 816 666 667$	24	116	116	24
$101.157 182 777 778$	24	116	116	24
$101.164 106 666 667$	24	116	116	24

Table 17. Table 17 : Success rate of outlier detection with real measurements

Case #	Dynamics	True positives	False negatives	False positives	True negatives
E,F	LF	13	5	0	0
E,F	MF	18	0	0	0
G,H	LF	10	5	0	3
G,H	MF	15	0	0	3

Table 18. Table 18 : True state and initial guesses for scenarios E to H

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
True	$- 21 551.184 664 630 2$	$14 404.866 452 074 8$	$- 1082.462 558 770 53$	$- 3.580 403 901 491 13$	$- 0.736 464 589 895 252$	$0.001 943 794 765 007 89$
Guess E	$- 21 547.072 614 259 1$	$14 404.656 186 620 9$	$- 1081.990 386 528$	$- 3.580 857 164 354 14$	$- 0.735 899 837 761 463$	$0.002 004 915 486 370 04$
Guess F	$- 21 551.886 616 748 6$	$14 407.416 553 400 4$	$- 1081.659 566 188 29$	$- 3.580 121 902 149 36$	$- 0.736 320 077 278 301$	$0.002 064 234 912 334 56$
Guess G	$- 21 547.072 614 259 1$	$14 404.656 186 620 9$	$- 1081.990 386 528$	$- 3.580 857 164 354 14$	$- 0.735 899 837 761 463$	$0.002 004 915 486 370 04$
Guess H	$- 21 551.886 616 748 6$	$14 407.416 553 400 4$	$- 1081.659 566 188 29$	$- 3.580 121 902 149 36$	$- 0.736 320 077 278 301$	$0.002 064 234 912 334 56$

Table 19. Table 19 : LS and LSAR estimation errors and ± 3 σ plus-or-minus 3 𝜎 \pm 3\sigma bounds for scenario E

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
LS Error	$9.239 450 350 005 2$	$3.305 740 148 502 32$	$0.116 199 529 730 007$	$0.001 753 121 538 532 74$	$0.000 306 294 804 796 094$	$0.000 136 508 035 645 701$
$\pm 3 σ$ bounds	$0.780 987 256 152 968$	$0.414 437 510 667 75$	$0.138 698 947 819 167$	$0.000 150 296 943 704 622$	$3.300 711 221 697 6 × 10^{- 05}$	$1.084 827 130 452 51 × 10^{- 05}$
LSAR Error	$0.602 359 107 738 526$	$2.323 684 898 290 35$	$0.774 951 866 708 533$	$0.000 183 186 873 129 161$	$7.607 963 996 845 18 × 10^{- 05}$	$0.000 171 030 451 399 098$
$\pm 3 σ$ bounds	$0.966 984 353 070 529$	$0.513 004 613 089 203$	$0.172 032 715 398 441$	$0.000 186 008 968 325 118$	$4.085 386 606 396 67 × 10^{- 05}$	$1.344 181 718 186 32 × 10^{- 05}$

Table 20. Table 20 : LS and LSAR estimation errors and ± 3 σ plus-or-minus 3 𝜎 \pm 3\sigma bounds for scenario G

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
LS Error	$464.561 341 965 927$	$289.620 372 040 554$	$31.661 500 496 537 3$	$0.093 766 813 232 237 1$	$0.020 533 983 440 204 7$	$0.001 565 373 075 141 34$
$\pm 3 σ$ bounds	$0.767 690 130 097 545$	$0.409 522 817 594 827$	$0.143 979 355 581 707$	$0.000 145 595 178 122 3$	$3.172 312 592 919 4 × 10^{- 05}$	$1.069 892 682 612 77 × 10^{- 05}$
LSAR Error	$1.255 601 626 539 59$	$2.693 725 983 598 98$	$0.819 188 049 874 167$	$0.000 288 182 607 168 274$	$2.516 861 659 601 47 × 10^{- 05}$	$0.000 180 207 964 448 427$
$\pm 3 σ$ bounds	$1.006 608 476 820 05$	$0.540 640 712 258 467$	$0.176 036 854 759 609$	$0.000 194 414 857 751 158$	$4.269 137 916 954 19 × 10^{- 05}$	$1.371 020 043 660 34 × 10^{- 05}$

Table 21. Table 21 : LS and LSAR estimation errors and ± 3 σ plus-or-minus 3 𝜎 \pm 3\sigma bounds for scenario H

	$x$ , $km$	$y$ , $km$	$z$ , $km$	$v_{x}$ , $km s^{- 1}$	$v_{y}$ , $km s^{- 1}$	$v_{z}$ , $km s^{- 1}$
LS Error	$3.296 045 131 536 92$	$0.326 729 321 715 63$	$0.505 709 718 794 199$	$0.000 549 673 851 591 107$	$4.288 232 319 844 32 × 10^{- 05}$	$0.000 156 495 686 685 633$
$\pm 3 σ$ bounds	$0.923 162 975 723 824$	$0.512 041 529 535 699$	$0.149 646 742 704 72$	$0.000 180 287 872 758 937$	$3.955 677 817 805 05 × 10^{- 05}$	$1.153 318 460 807 × 10^{- 05}$
LSAR Error	$1.151 814 072 398 7$	$2.635 204 221 686 35$	$0.815 017 315 887 631$	$0.000 268 448 717 855 625$	$2.938 106 897 494 42 × 10^{- 05}$	$0.000 179 925 964 703 669$
$\pm 3 σ$ bounds	$1.144 552 955 839 21$	$0.634 817 164 172 961$	$0.185 595 487 082 44$	$0.000 223 489 960 097 547$	$4.903 226 099 498 16 × 10^{- 05}$	$1.429 585 524 512 57 × 10^{- 05}$

Table 22. Table 22 : Runtime and number of iterations with real measurements

Case #	Pruning	LS	LSAR	LS	LSAR
	$t$ , $s$			# iterations
E	-	$510.904$	$200.603$	3	3
F	$49.653$	$453.881$	$198.409$	3	3
G	-	$537.406$	$192.47$	5	3
H	$14.897$	$464.756$	$140.941$	5	3

Equations120

\overset{x}{˙} = f (x, t) + B (t) w (t)

\overset{x}{˙} = f (x, t) + B (t) w (t)

E [w (t)]

E [w (t)]

E [w (t) w^{T} (t)]

{\overset{x}{˙} (t) \dot{P} (t) = f (x, t) = A (t) P (t) + P (t) A^{T} (t) + B (t) Q (t) B^{T} (t)

{\overset{x}{˙} (t) \dot{P} (t) = f (x, t) = A (t) P (t) + P (t) A^{T} (t) + B (t) Q (t) B^{T} (t)

r = r_{o b s} + ρ \overset{ρ}{^}

r = r_{o b s} + ρ \overset{ρ}{^}

\overset{ρ}{^} = cos α cos δ sin α cos δ sin δ

\overset{ρ}{^} = cos α cos δ sin α cos δ sin δ

⎩ ⎨ ⎧ ρ = ρ_{x}^{2} + ρ_{y}^{2} + ρ_{z}^{2} α = tan^{- 1} ρ_{y} / ρ_{x} δ = sin^{- 1} ρ_{z} / ρ

⎩ ⎨ ⎧ ρ = ρ_{x}^{2} + ρ_{y}^{2} + ρ_{z}^{2} α = tan^{- 1} ρ_{y} / ρ_{x} δ = sin^{- 1} ρ_{z} / ρ

= α_{i} + c σ_{α_{i}} \cdot δ α_{i}

= α_{i} + c σ_{α_{i}} \cdot δ α_{i}

[δ_{i}]

f \approx [f] = T_{f} (δ x)

f \approx [f] = T_{f} (δ x)

b_{l}

b_{l}

b_{u}

p (δ x) \in [a_{0} \pm k = 1 \sum v \absolutevalue a_{k}^{o}]

p (δ x) \in [a_{0} \pm k = 1 \sum v \absolutevalue a_{k}^{o}]

[x] = \overset{x}{ˉ} + β ⊙ δ x

[x] = \overset{x}{ˉ} + β ⊙ δ x

[y] = f ([x]) = T_{y} (δ x)

[y] = f ([x]) = T_{y} (δ x)

[J] = \frac{\partial f}{\partial x} = \overline{J} + δ J

[J] = \frac{\partial f}{\partial x} = \overline{J} + δ J

b_{ij} = k = 1 \sum n \absolutevalue a_{ij, k}

b_{ij} = k = 1 \sum n \absolutevalue a_{ij, k}

ν = \frac{\sum _{i = 1}^{m} \sum _{j = 1}^{n} ( \sum _{k = 1}^{n} \absolutevalue a _{ij, k} ) ^{2}}{\sum _{i = 1}^{m} \sum _{j = 1}^{n} J _{ij}^{2}}

ν = \frac{\sum _{i = 1}^{m} \sum _{j = 1}^{n} ( \sum _{k = 1}^{n} \absolutevalue a _{ij, k} ) ^{2}}{\sum _{i = 1}^{m} \sum _{j = 1}^{n} J _{ij}^{2}}

[J_{d}] = [J] \circ δ x_{d}

[J_{d}] = [J] \circ δ x_{d}

d_{s} = d arg max {ν_{d}}

d_{s} = d arg max {ν_{d}}

[x^{(j)}] = [x^{(0)}] \circ {δ x_{1}, \dots, \frac{2}{3} \cdot (j - 2) + \frac{1}{3} \cdot δ x_{d_{s}}, \dots, δ x_{v}}^{T}

[x^{(j)}] = [x^{(0)}] \circ {δ x_{1}, \dots, \frac{2}{3} \cdot (j - 2) + \frac{1}{3} \cdot δ x_{d_{s}}, \dots, δ x_{v}}^{T}

[z^{(i)}] = [z^{(2)}] \circ {δ x_{1}, \dots, 3 δ x_{d_{s}}, \dots, δ x_{v}}^{T}

[z^{(i)}] = [z^{(2)}] \circ {δ x_{1}, \dots, 3 δ x_{d_{s}}, \dots, δ x_{v}}^{T}

\overset{y}{^}_{k} = h_{k} (\overset{x}{^}_{k}) k \in [1, m]

\overset{y}{^}_{k} = h_{k} (\overset{x}{^}_{k}) k \in [1, m]

h_{k} (x_{k})

h_{k} (x_{k})

= h_{k} (\overset{x}{^}_{k}) + \tilde{H}_{k} Φ (t_{k}, t_{0}) δ x_{0}

= h_{k} (\overset{x}{^}_{k}) + H_{k} δ x_{0}

\tilde{H}_{k} = \frac{\partial h _{k}}{\partial x _{k}}_{\overset{x}{^}_{k}}

\tilde{H}_{k} = \frac{\partial h _{k}}{\partial x _{k}}_{\overset{x}{^}_{k}}

H_{k} = \tilde{H}_{k} Φ (t_{k}, t_{0}) = \frac{\partial h _{k}}{\partial x _{k}} \frac{\partial x _{k}}{\partial x _{0}}_{\overset{x}{^}_{0}}

e

e

= y - h (\overset{x}{^}) - H δ x_{0}

= Δ y - H δ x_{0}

y

y

\overset{y}{^}

h (\overset{x}{^})

H

\overset{x}{^}_{0}^{j} = \overset{x}{^}_{0}^{j - 1} + δ x_{0}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsTarget Tracking and Data Fusion in Sensor Networks · Fault Detection and Control Systems · Inertial Sensor and Navigation

Full text

\setabbreviationstyle

[acronym]long-short

A Multifidelity Approach to Robust Orbit Determination

Alberto Fossà

Department of Aerospace Vehicles Design and Control, Institut Supérieur de l’Aéronautique et de l’Espace, 10 Avenue Edouard Belin, Toulouse, 31055, France

Corresponding author, [email protected]

Roberto Armellin

Te Pūnaha Ātea - Space Institute, The University of Auckland, 20 Symonds Street, Auckland, 1010, New Zealand

Emmanuel Delande

Space Situational Awareness Office, Centre National d’Études Spatiales, 18 Avenue Edouard Belin, Toulouse, 31400, France

Matteo Losacco

Department of Aerospace Vehicles Design and Control, Institut Supérieur de l’Aéronautique et de l’Espace, 10 Avenue Edouard Belin, Toulouse, 31055, France

Francesco Sanfedino

Department of Aerospace Vehicles Design and Control, Institut Supérieur de l’Aéronautique et de l’Espace, 10 Avenue Edouard Belin, Toulouse, 31055, France

Abstract

This paper presents an algorithm for the preprocessing of observation data aimed at improving the robustness of orbit determination tools. Two objectives are fulfilled: obtain a refined solution to the initial orbit determination problem and detect possible outliers in the processed measurements. The uncertainty on the initial estimate is propagated forward in time and progressively reduced by exploiting sensor data available in said propagation window. Differential algebra techniques and a novel automatic domain splitting algorithm for second-order Taylor expansions are used to efficiently propagate uncertainties over time. A multifidelity approach is employed to minimize the computational effort while retaining the accuracy of the propagated estimate. At each observation epoch, a polynomial map is obtained by projecting the propagated states onto the observable space. Domains that do no overlap with the actual measurement are pruned thus reducing the uncertainty to be further propagated. Measurement outliers are also detected in this step. The refined estimate and retained observations are then used to improve the robustness of batch orbit determination tools. The effectiveness of the algorithm is demonstrated for a geostationary transfer orbit object using synthetic and real observation data from the TAROT network.

Keywords: uncertainty propagation, multifidelity methods, robust orbit determination

1 Introduction

Preserving the realism of the probabilistic representation of space objects (SOs) states over time is paramount in space situational awareness (SSA) applications such as orbit determination (OD) and SOs catalogs build-up and maintenance. To accomplish these tasks, two main tools are required: an accurate and efficient uncertainty propagation (UP) method and an effective strategy to incorporate newly available information in the predicted estimate. This work contributes to these needs by proposing a pruning algorithm for the refinement of an initial solution, typically the result of an initial orbit determination (IOD) process, and the detection of measurements outliers among the input set. The algorithm is composed of four main steps which are carried out sequentially for each available measurements: propagation of the current state estimate, projection onto the observable space, pruning of the uncertainty based on the newly available information, and merging to reduce the information to be further propagated.

Regarding the UP problem in orbital dynamics, this has been historically tackled with either linearized methods or large scale Monte Carlo (MC) simulations. Several methods for nonlinear UP have since been proposed and are conveniently classified as intrusive or non-intrusive, with the first requiring full knowledge of the underlying dynamics while the seconds treating the dynamical model as a black box. The first category includes State Transition Tensors (STTs) [1] and differential algebra (DA) [2, 3] while the second gathers the aforementioned MC as well as unscented transform (UT) [4], conjugate unscented transform (CUT) [5] and polynomial chaos expansions (PCEs) [6]. A comprehensive survey of these methods can be found in [7]. In this work, DA is used for UP in conjunction with an automatic domain splitting (ADS) algorithm [8] to control the error on the final solution manifold. Multifidelity methods are exploited to drastically reduce the computational effort of the overall algorithm while retaining a similar accuracy on the final solution. Their effectiveness for orbit UP was already demonstrated in [9, 10] and corroborated in this work. The DA-based UP is carried out on a low-fidelity dynamical model and only the centers of the resulting manifold are propagated in high-fidelity to compute a posteriori correction to the low-fidelity solution.

The projection onto the observables space follows, also wrapped by the ADS algorithm to retain the accuracy of the resulting polynomial map, or manifold, that represents the uncertainty on the predicted measurements by patching a set of truncated Taylor expansions. Polynomial bounding techniques are then employed to estimate the range of variation of both predicted and actual measurements. These bounds are then used to discard the domains that do not match the real observation and reduce the size of the manifold of states to be further propagated. Possible measurements outliers are also identified in the observations that have no intersection with the projected set. The output is a refined solution to the IOD problem that is used as improved guess to initialize a batch OD algorithm. Concerning the last, if least squares (LS) techniques are the workhorse of current operational algorithms, estimators that minimize the $L_{1}$ norm of the residuals vector were shown to improve robustness in the presence of measurements outliers [11]. These () estimators were also proven to be faster than LS regression on large data sets [12] and successfully applied to the problem of batch OD [13, 14]. The least sum of absolute residuals (LSAR) algorithm was then proposed in [15] as an alternative technique to compute the minimum $L_{1}$ norm estimates. This algorithm recasts the problem into a sequence of linear programming problems (LPPs) which can be then solved with highly efficient simplex or interior point methods [16]. Both solutions are compared in this paper highlighting the benefits provided by the pruning algorithm over a direct processing of all available information. A complete solution to the batch OD problem is thus presented, with DA being the common denominator underlying each proposed tool. This approach is demonstrated of being robust to measurements outliers thanks to an effective preprocessing of the available observations and the use of the LSAR algorithm to compute the final estimate.

The paper is organized as follows. Section 2 sets the background material for the paper. This includes the modeling of the problem in Section 2.1, the DA framework in Section 2.2 and the OD tools in Section 2.3. The main contribution of the current work, namely the multifidelity pruning scheme, is then detailed in Section 3 after the DA-based IOD algorithm is briefly outlined. Numerical applications of the proposed tool are presented in Section 4 before drawing the conclusions in Section 5.

2 Mathematical background

This section describes the building blocks this work is built on. These include the modeling of the process noise and the measurement process, differential algebra techniques applied to uncertainty propagation, the automatic domain splitting and merging algorithms used to control the accuracy of Taylor expansions, and the batch orbit determination tools.

2.1 Modeling of the problem

The numerical algorithm used to model the process noise is introduced hereafter, followed by a description of the measurement process and sensor noise model for optical telescopes.

2.1.1 State noise compensation

Taking into account the uncertainty due to neglected or mismodeled accelerations is paramount to achieve a precise estimate of the propagated orbit. Stochastic accelerations are thus included in the developed UP framework and their effects on the state covariance estimated using the State Noise Compensation (SNC) algorithm [17].

For continuous systems, the dynamics under the influence of process noise is described by

[TABLE]

with ${\bf\it f}:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$ the deterministic dynamics, ${\bf\it w}(t)$ the $m$ -dimensional process noise, and ${\bf\it B}(t)$ the $n\times m$ process noise mapping matrix that maps the process noise into state derivatives. Under the assumption of uncorrelated white noise, the first two statistical moments of ${\bf\it w}(t)$ are given by

[TABLE]

with $\delta(t-\tau)$ the Dirac delta and ${\bf\it Q}(t)$ the function of ${\bf\it w}(t)$ [18]. An estimate of the true state is then maintained through a mean state ${\bf\it x}(t)$ and a covariance matrix ${\bf\it P}(t)$ that are computed by integrating the system of first-order ordinary differential equations (ODEs)

[TABLE]

subject to initial conditions ${\bf\it x}(t_{0})={\bf\it x}_{0}$ and ${\bf\it P}(t_{0})={\bf\it 0}$ .

2.1.2 Measurement process and sensor noise

Even though the pruning algorithm presented in Section 3.2 is agnostic with respect to the type of measurements, simulations in Section 4 uses angular measurements from optical sensors. The corresponding measurement model is thus briefly introduced hereafter.

Consider an optical telescope identified by its geodetic coordinates $(\phi,\lambda,h)$ with $\phi\in[-\pi/2,\pi/2]$ geodetic latitude, $\lambda\in[-\pi,\pi]$ geodetic longitude and $h\in\mathbb{R}^{+}$ geodetic height. Denote with ${\bf\it x}=[{\bf\it r}^{T}\ {\bf\it v}^{T}]^{T}$ the inertial state vector of the tracked SO and with ${\bf\it r}_{obs}$ the inertial position of the telescope at the same epoch. The two positions ${\bf\it r},{\bf\it r}_{obs}$ are related by

[TABLE]

with $\rho$ topocentric range and $\hat{{\bf\it\rho}}$ line of sight (LOS) unit vector. The last is computed from the topocentric right ascension and declination $(\alpha,\delta)$ as

[TABLE]

If the LOS vector $\rho\hat{{\bf\it\rho}}=[\rho_{x}\ \rho_{y}\ \rho_{z}]^{T}$ is known instead, the topocentric range, right ascension and declination are given by

[TABLE]

Sensor noise is modeled by assuming uncorrelated angular measurements $(\alpha_{i},\delta_{i})$ corrupted by Gaussian white noise with standard deviations $(\sigma_{\alpha_{i}},\sigma_{\delta_{i}})$ for each measurement epoch $t_{i}$ . These uncertainties are represented in the DA framework as

[TABLE]

where $c>0$ is the $z$ -score that quantifies how many standard deviations lie within the polynomial bounds. This parameter is typically set equal to $c=3$ such that the true measurement is found within the domain of the polynomial with a $\approx 99.73\%$ ( $3\sigma$ ) confidence.

2.2 Differential algebra

Differential algebra is a computing technique based on the idea that it is possible to extract more information from a function $f$ than its mere value $y=f(x)$ at $x$ . It replaces the algebra of floating point numbers with a new algebra of Taylor polynomials to compute the given order Taylor expansion of any function $f$ in $v$ variables that is $\mathcal{C}^{n+1}$ in the domain of interest $[-1,1]^{v}$ [2]. The notation used is

[TABLE]

where $\delta{\bf\it x}=\{\delta x_{1},\ldots,\delta x_{v}\}^{T}$ are the $v$ independent DA variables.

As for the floating point algebra, the four arithmetic operations ( $+,-,\ast,\div$ ), elementary functions (e.g. trigonometric functions, powers, exponential and logarithm), derivation, integration, function composition and inversion are well defined in DA. The latter technique is particularly suitable for the solution of implicit equations and the computation of the Taylor expansion of the flow of differential equations in terms of their initial conditions.

2.2.1 Polynomial bounding

Consider the DA of order $n$ in $v$ variables, and a multivariate polynomial $p(\delta{\bf\it x})$ with coefficients $a_{k}\in\mathbb{R}$ . Define as odd the coefficients for which at least one DA variable appears with an odd exponent (e.g. $\delta x_{1}$ or $\delta x_{1}^{2}\delta x_{2}$ ), and denote them with the superscript “ $o$ ”. Likewise, define as even coefficients those for which all variables appear with even exponents (e.g $\delta x_{1}^{2}$ or $\delta x_{1}^{2}\delta x_{2}^{2}$ ), and denote them with the superscript “ $e$ ”. Since $\delta{\bf\it x}\in[-1,1]^{v}$ , an estimation of the range of variation, or bounds, of $p(\delta{\bf\it x})$ can be obtained as

[TABLE]

with $b_{l},b_{u}$ lower and upper bounds such that $p(\delta{\bf\it x})\in[b_{l},b_{u}]$ , $a_{0}$ constant part and $n_{o},n_{e}$ number of odd and even coefficients respectively. For first-order polynomials the bounds computed with Eq. 9 are rigorous in the sense that they coincide with the extrema of $p(\delta{\bf\it x})$ when evaluated in $[-1,1]^{v}$ . In this case all even coefficients are identically zero and Eq. 9 simplifies to

[TABLE]

2.2.2 Uncertainty propagation

Given the nonlinear transformation ${\bf\it f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ and two multivariate random variables ${\bf\it X},{\bf\it Y}$ with probability density function (pdf) $p_{{\bf\it X}},p_{{\bf\it Y}}$ respectively, the UP problem requires the estimation of $p_{{\bf\it Y}}$ from the only knowledge of $p_{{\bf\it X}}$ and ${\bf\it f}$ . In DA framework, ${\bf\it X}$ is represented as a vector of first-order multivariate polynomials in the form

[TABLE]

with $\bar{{\bf\it x}}=\mathbb{E}[{\bf\it X}]$ the expected value of ${\bf\it X}$ , ${\bf\it\beta}$ the semi-amplitude of the range of variation of $[{\bf\it x}]$ , $\delta{\bf\it x}$ the first-order variations of ${\bf\it X}$ around $\bar{{\bf\it x}}$ and $\odot$ the Hadamard (or element-wise) product. The Taylor expansion of ${\bf\it Y}$ in terms of $\delta{\bf\it x}$ is then obtained evaluating ${\bf\it f}$ in DA framework as

[TABLE]

In the context of UP, different methods were proposed to efficiently compute the statistical properties of ${\bf\it Y}$ from $\mathcal{T}_{{\bf\it y}}(\delta{\bf\it x})$ thus avoiding the needs for expensive Monte Carlo (MC) simulations [3].

2.2.3 Low-order automatic domain splitting

Taylor polynomials are only a local approximation of the function ${\bf\it f}$ around the expansion point, and a single polynomial might fail to accurately represent the uncertainty for $\absolutevalue{\delta{\bf\it x}}\gg 0$ . To overcome this limitation, automatic domain splitting (ADS) was firstly proposed by [8] to monitor the accuracy of the expansion at the edges of its domain and recursively split the initial polynomial into smaller subdomains in which the required accuracy is satisfied.

Instead, a novel splitting technique specific for second-order Taylor polynomials is used in this work [10, 19]. The algorithm, named low-order automatic domain splitting (LOADS), defines a nonlinearity index (NLI) computed from the Jacobian of the transformation $\partial{\bf\it f}/\partial{\bf\it x}$ that was inspired by the nonlinearity measure introduced by [20]. In this paper, ${\bf\it f}$ is referred as the target function of the LOADS scheme.

Given an initial domain $[{\bf\it x}^{(0)}]$ , the algorithm starts by evaluating ${\bf\it f}$ in $[{\bf\it x}^{(0)}]$ and computing the NLI of the output $[{\bf\it y}^{(0)}]$ as follows.

The Jacobian of the transformation is firstly computed as the sum of its constant part $\overline{{\bf\it J}}$ and first-order expansion $\delta{\bf\it J}(\delta{\bf\it x})$ as

[TABLE]

since $\delta{\bf\it J}\equiv{\bf\it 0}$ for a linear transformation, a measure of nonlinearity can be obtained from the magnitude of $\delta{\bf\it J}$ . For each entry $\delta J_{ij}$ , rigorous bounds are firstly obtained from Eq. 9 as

[TABLE]

where $a_{ij,k}$ are the first-order coefficients of $J_{ij}$ , i.e. $\delta J_{ij}=\sum_{k=1}^{n}a_{ij,k}\delta x_{k}$ and thus $\delta J_{ij}\in[-b_{ij},b_{ij}]$ . The nonlinearity index is then defined as

[TABLE]

which is the ratio between the Frobenius norm of the matrix of the bounds and the Frobenius norm of the constant part of the Jacobian $\overline{{\bf\it J}}$ .

The computed NLI is then compared with an imposed error threshold $\varepsilon$ and a decision is made. If $\nu\leq\varepsilon$ the algorithm ends and $[{\bf\it y}^{(0)}]$ is retained. On the contrary, if $\nu>\varepsilon$ , the output is rejected and $[{\bf\it x}^{(0)}]$ is split into three polynomials $[{\bf\it x}^{(j)}]$ with each one covering $1/3$ of the initial domain.

To perform the split, a second index is firstly computed to identify the direction that contributes the most to the nonlinearity in $[{\bf\it y}]$ . For each independent DA variable $\delta x_{d},\ d\in[1,v]$ , a directional Jacobian is built as the function composition of $[{\bf\it J}]$ and $\delta{\bf\it x}_{d}$ as

[TABLE]

where $\delta{\bf\it x}_{d}$ is a vector of zeros except for the $d^{th}$ component set equal to $\delta x_{d}$ , i.e. $\delta{\bf\it x}_{d}=\{0,\ldots,\delta x_{d},\ldots,0\}^{T}$ , and $\circ$ denotes function composition. A directional NLI $\nu_{d}$ is then computed from $[{\bf\it J}_{d}]$ as in Eq. 15 and the splitting direction $d_{s}$ selected for the maximum $\nu_{d}$

[TABLE]

Knowing $d_{s}$ , the three subdomains $[{\bf\it x}^{(j)}]$ for $j\in[1,3]$ are computed from $[{\bf\it x}^{(0)}]$ as

[TABLE]

The target function is then evaluated on the set $\{[{\bf\it x}^{(j)}]\}$ and the procedure iterated until the accuracy is satisfied in each subdomain. The result is a set of polynomials for both the initial domain and its image through ${\bf\it f}$ . These sets are named manifolds and denoted as $M_{{\bf\it x}}=\{[{\bf\it x}^{(i)}]\}$ and $M_{{\bf\it y}}=\{[{\bf\it y}^{(i)}]\}$ respectively where $i\in[1,N]$ and $N$ is the total number of domains after splitting.

2.2.4 Domain merging

Consider a second function ${\bf\it g}:\mathbb{R}^{m}\to\mathbb{R}^{p}$ evaluated on the manifold $M_{{\bf\it y}}$ to obtain $M_{{\bf\it z}}$ . Nonlinearities in $[{\bf\it z}^{(i)}]$ can be smaller than that of $[{\bf\it y}^{(i)}]$ and one or more \glsxtrshortplnli computed from $\partial{\bf\it g}/\partial{\bf\it y}$ can fall far below the splitting threshold, i.e. $\nu^{(i)}\ll\varepsilon$ . In this case is thus conceivable to develop an algorithm to recombine the polynomials $[{\bf\it z}^{(i)}]$ and reduce the size of $M_{{\bf\it z}}$ . The merging scheme used in this work stems from this observation and requires the knowledge of the entire splitting history of each polynomial in $M_{{\bf\it z}}$ to attempt a recombination of the different domains [19]. This information comprises all directions $d_{s}$ and shifts $j$ along which $[{\bf\it z}^{(i)}]$ was split and allows to identify domains that generated from the same parent and thus potentially mergeable.

Given the manifold $M_{{\bf\it z}}$ , the merging scheme proceeds as follows. Firstly, the domains are grouped by depth of split, i.e. by the number of times they have been split. Secondly, starting from the sub-manifold with the greatest depth, the histories of its elements are compared to identify all triplets sharing a common parent. For each triplet, the domain of the parent is recovered as

[TABLE]

with $d_{s}$ last splitting direction and $[{\bf\it z}^{(2)}]$ central polynomial. A decision on the potential merge is then taken after computing the NLI of $[{\bf\it z}^{(i)}]$ : if $\nu\leq\varepsilon$ , $[{\bf\it z}^{(i)}]$ is retained and added to the sub-manifold with a shallower depth. If $\nu>\varepsilon$ , $[{\bf\it z}^{(i)}]$ is discarded and the three domains of the original triplet added to the final manifold $M^{\prime}_{{\bf\it z}}$ . This procedure is repeated until all sub-manifolds have been emptied. At each stage, if an incomplete triplet is found, no merge is attempted and its elements are directly added to $M^{\prime}_{{\bf\it z}}$ .

2.3 Orbit determination

The problem of orbit determination consists in estimating the state of an orbiting object given the measurements provided by one or more sensors. Among the different approaches to the OD problem, this paper exploits the so-called batch algorithms introduced in the followings.

2.3.1 Batch orbit determination algorithms

Solving the batch OD problem requires to find an estimate of the orbit state $\hat{{\bf\it x}}_{0}=\hat{{\bf\it x}}(t_{0})$ , at some reference epoch $t_{0}$ , that minimizes the residuals between the real observations $({\bf\it y}_{k})_{k=1}^{m}$ collected at epochs $(t_{k})_{k=1}^{m}$ and belonging to some observation space $\mathbf{Y}\subseteq\mathbb{R}^{p}$ , and the expected (or theoretical) measurements $(\hat{{\bf\it y}}_{k})_{k=1}^{m}$ . The latter are computed by propagating $\hat{{\bf\it x}}_{0}$ to $t_{k}$ and projecting the state $\hat{{\bf\it x}}_{k}$ onto the observable space, i.e.,

[TABLE]

where ${\bf\it h}_{k}$ is the observation function at epoch $t_{k}$ . Assuming that the estimated state $\hat{{\bf\it x}}_{k}$ is close to the true object’s state ${\bf\it x}_{k}$ at epoch $t_{k}$ , ${\bf\it h}_{k}$ can be linearized using a first-order Taylor expansion as

[TABLE]

with ${\bf\it\Phi}(t_{k},t_{0})$ the state transition matrix (STM) from $t_{0}$ to $t_{k}$ , $\tilde{{\bf\it H}}_{k}$ the partials of ${\bf\it h}_{k}$ with respect to ${\bf\it x}_{k}$ and evaluated at $\hat{{\bf\it x}}_{k}$ , and ${\bf\it H}_{k}$ the partials of ${\bf\it h}_{k}$ with respect to ${\bf\it x}_{0}$ and evaluated at $\hat{{\bf\it x}}_{0}$ . The last two matrices are defined as

[TABLE]

The residuals are then readily computed as

[TABLE]

where the total observations vectors and Jacobian are defined as

[TABLE]

The three vectors ${\bf\it y},\hat{{\bf\it y}}$ and ${\bf\it h}(\hat{{\bf\it x}})$ have size $mp\times 1$ while the matrix ${\bf\it H}$ has size $mp\times n$ , where $n$ is the dimension of $\hat{{\bf\it x}}_{k}$ .

Given an initial guess for $\hat{{\bf\it x}}_{0}$ denoted as $\hat{{\bf\it x}}_{0}^{0}$ , the objective of the OD algorithm is to iteratively solve for $\delta{\bf\it x}_{0}$ and compute an updated estimate $\hat{{\bf\it x}}_{0}^{j}$ as

[TABLE]

with $j$ iteration number. The recursive procedure terminates as soon as one of the following criteria is satisfied

•

$\norm{\hat{{\bf\it y}}^{j}-\hat{{\bf\it y}}^{j-1}}_{2}<\varepsilon_{res}$ , i.e. the residuals at two subsequent iterations differ by less than the residual tolerance $\varepsilon_{res}$ .

•

$\absolutevalue{J^{j}-J^{j-1}}<\varepsilon_{opt}$ , i.e. the cost functions at two subsequent iterations differ by less than the optimality tolerance $\varepsilon_{opt}$ .111the cost function $J$ is defined in Sections 2.3.2 and 2.3.3.

•

$\norm{\delta\hat{{\bf\it x}}_{0}}_{2}<\varepsilon_{step}$ , i.e. the state update is smaller than the step tolerance $\varepsilon_{step}$ .

•

$j=j_{max}$ , i.e. a maximum number of iterations has been reached.

2.3.2 Least squares

The LS solution to the OD problem solves for $\delta{\bf\it x}_{0}$ that minimizes the residuals in the LS sense. The cost function is given by

[TABLE]

with ${\bf\it W}=\text{diag}(1/\sigma_{i}^{2})$ the diagonal weight matrix where $\sigma_{i}$ is the sensor noise standard deviation corresponding to the $i^{th}$ observation. The state update that minimizes Eq. 26 is obtained using the Levenberg–Marquardt algorithm [21] (damped LS) as

[TABLE]

with $\lambda$ the damping parameter and $I$ the identity matrix. The covariance matrix of $\hat{{\bf\it x}}_{0}$ is finally estimated as

[TABLE]

2.3.3 Minimum $L_{1}$ norm estimates

The minimum $L_{1}$ norm solution to the OD problem was proven to be a robust alternative to LS optimization capable of rejecting outlier measurements in the presence of redundant observations [11, 12, 13, 14, 15]. In this work, the least sum of absolute residuals (LSAR) algorithm proposed in [15] is used to solve the OD problem in a $L_{1}$ sense. The cost function to be minimized is

[TABLE]

with $e_{i}$ the components of ${\bf\it e}$ and $\tilde{w}_{i}$ the weights defined as $\tilde{w}_{i}=1/(1.24\sigma_{i})$ [22]. The problem is recast as a linear programming problem (LPP) in the form

[TABLE]

with

[TABLE]

and ${\bf\it s}$ slack variables for the constraint vector ${\bf\it e}$ such that

[TABLE]

Throughout this work, the solver MOSEK222https://www.mosek.com with its interior point algorithm is used to efficiently solve Eq. 30. The covariance matrix of $\hat{{\bf\it x}}_{0}$ is then estimated as

[TABLE]

with $\tilde{{\bf\it W}}=\text{diag}(\tilde{w}_{i}^{2})={\bf\it W}/1.24^{2}$ .

3 Robust orbit determination

The full batch orbit determination chain is presented in this section. Given the space object (SO) of interest, an estimate of its orbit state at epoch $t_{0}$ is firstly obtained with the DA-based initial orbit determination (IOD) routine described in Section 3.1. This estimate is then propagated forward in time using multifidelity techniques to minimize the computational effort and its uncertainty progressively reduced with a custom pruning scheme as soon as new measurements become available. Potential outliers in the input measurements are also identified and discarded in this step. A batch OD algorithm is finally run on the pruned state and retained observations to compute the final estimate. Both LS and LSAR techniques are tested and their performance discussed in Section 4. The overall workflow is summarized in Fig. 1. This paper assumes observations from optical telescopes, but each tool can be easily adapted to other types of measurements. For the IOD algorithm in Section 3.1 it is already the case with a recently submitted paper that extends the algorithm to radar measurements [23].

3.1 Initial orbit determination

Consider a set of $N$ tuples of angular measurements as provided by a ground based optical sensor while observing an unknown SO

[TABLE]

with $\alpha_{i},\delta_{i}$ topocentric right ascension and declination of the SO at epoch $t_{i}$ and $\sigma_{\alpha_{i}},\sigma_{\delta_{i}}$ associated standard deviations of the sensor noise, assumed as uncorrelated Gaussian white noise as described in Section 2.1.2.

The goal of the IOD is to estimate the state of the object from the available set of measurements. The problem is well-known in literature and typically solved by processing the nominal angular measurements only [24, 25]. Instead, the approach here exploited makes use of DA to obtain an analytic map between the confidence intervals of the processed measurements and the uncertainty of the computed solution. It was firstly developed under the assumption of Keplerian dynamics by [26] for optical measurements and [27] for radar sensors, and recently extended in [23] to accommodate any dynamical model for the target object. The method sets up an iterative procedure to estimate the topocentric ranges of the object at the epoch of the first, middle, and last available measurements, namely $\rho_{k}=\{\rho_{1},\rho_{2},\rho_{3}\}$ . The initial guess for these ranges is provided by the classical Gauss solution [25]. Two Lambert’s problems [28] are then solved between $t_{1},t_{2}$ and $t_{2},t_{3}$ respectively. A nonlinear solver is then used to update $\rho_{j}$ such that the velocity vectors at the intermediate epoch coincide. An estimate of the object state is then directly available from this solution.

A polynomial expansion of the state with respect to the observables can be obtained by initializing the processed measurements as Taylor variables and performing all operations in the DA framework. Under the assumption of uncorrelated measurements corrupted by Gaussian white noise, these are written as

[TABLE]

where $c>0$ is the $z$ -score introduced in Section 2.1.2 and

[TABLE]

The state vectors at $t_{j}$ can then be written as

[TABLE]

The solution in Eq. 37 is obtained under the assumption of unperturbed Keplerian motion (since based on the recursive solution of two Lambert problems) and might be inaccurate for long observation windows in which deviations from the nominal two-body trajectory cannot be neglected. A DA-based shooting scheme is thus employed to account for the effects of the $J_{2}$ term of the Earth’s gravitational potential. Given the needs for a fast propagation method at this stage, the analytical formulation of the $J_{2}$ -perturbed dynamics proposed in [29] is employed. Starting from Eq. 37, an updated estimate is then obtained as

[TABLE]

To maintain an accurate polynomial representation of the state in the presence of large measurements errors, the IOD algorithm is wrapped within the LOADS framework introduced in Section 2.2.3. A manifold of polynomials $M_{{\bf\it x}_{j}^{J_{2}}}$ is thus automatically generated while computing the expansion of the IOD solution with respect to $({\bf\it\alpha},{\bf\it\delta})$ such that the nonlinear transformation is accurately described across its entire domain. These manifolds are defined as

[TABLE]

where $N$ is the number of domains generated by the LOADS algorithm and $j$ denotes the observation epoch.

3.2 Multifidelity uncertainty propagation

Multifidelity methods are introduced to improve the computational efficiency of the prediction step while guaranteeing the required accuracy for the final estimate. In particular, a bifidelity technique based on two distinct dynamical models is employed.

The low-fidelity model is the Simplified General Perturbations (SGP4) model used by NORAD to produce daily Two-Line Elements (TLEs) of all tracked SOs [30]. Given its analytical formulation, SGP4 is easily embedded in the DA-LOADS framework to efficiently propagate orbit states and associated uncertainties as Taylor polynomials. A high-fidelity correction is yet needed to retain the accuracy of the final solution. A numerical propagator is thus used to integrate the high-fidelity orbital dynamics as well as the effects of process noise via Eq. 3. In this work, ${\bf\it f}({\bf\it x},t)$ are the Gauss variational equations which describe the motion of SOs in the perturbed Keplerian dynamics [31]. A representation of the orbit state at epoch $t$ is then maintained through a manifold of Taylor polynomials $M_{MF}(t)=\{[{\bf\it x}_{MF}^{(i)}(t)]\}$ and a manifold of covariance matrices $M_{PN}(t)=\{{\bf\it P}_{PN}^{(i)}(t)\}$ characterized by a one-to-one correspondence between the elements in $M_{MF}(t)$ and the elements in $M_{PN}(t)$ . The polynomials $[{\bf\it x}_{MF}^{(i)}(t)]$ represent the initial uncertainty on the IOD solution propagated at time $t$ under a deterministic dynamical model and are centered on the high-fidelity reference trajectories solution of the first of Eq. 3. The covariance matrices ${\bf\it P}_{PN}^{(i)}(t)$ represent instead the contribution of process noise effects and are obtained by numerically integrating the second of Eq. 3. At initial time $t_{0}$ the first manifold $M_{MF}(t_{0})$ is thus given by Eq. 38 for $j=1$ while all covariance matrices are initialized to zero, i.e. ${\bf\it P}_{PN}^{(i)}(t_{0})={\bf\it 0}\ \forall i$ . Note that, at any time $t$ , the polynomials $[{\bf\it x}_{MF}^{(i)}(t)]$ are function of the six independent DA variables $(\delta{\bf\it\alpha},\delta{\bf\it\delta})$ , i.e.

[TABLE]

Consider the propagation time span $[t_{k-1},t_{k}]$ between two subsequent observations, and denote with $M_{MF}(t_{k-1})$ and $M_{PN}(t_{k-1})$ the manifolds of Taylor polynomials and covariance matrices at time $t_{k-1}$ respectively. The solution manifolds $M_{MF}(t_{k}),M_{PN}(t_{k})$ at time $t_{k}$ are then computed as follows. Firstly, $M_{MF}(t_{k-1})$ is processed with the LOADS scheme setting SGP4 as target function ${\bf\it f}$ . The results is a low-fidelity solution manifold $M_{LF}(t_{k})$ whose size is generally different from that of $M_{MF}(t_{k-1})$ . In the second step, the centers of the polynomials $[{\bf\it x}^{(i)}(t_{k-1})]$ that maps to $[{\bf\it x}^{(i)}(t_{k})]$ in $M_{LF}(t_{k})$ are propagated point-wise in high-fidelity using Eq. 3 to obtain the reference solutions $\{{\bf\it x}^{(i)}_{HF}(t_{k})\}$ and covariance matrices $\{{\bf\it P}^{(i)}_{PN}(t_{k})\}$ for each polynomial in $M_{LF}(t_{k})$ . Note that Eq. 3 is now subject to initial conditions $\bar{{\bf\it x}}^{(i)}_{MF}(t_{k-1}),{\bf\it P}^{(i)}_{PN}(t_{k-1})$ , while ${\bf\it P}^{(i)}_{PN}(t_{k})$ is the result of the effects of process noise and do not include the contribution of the initial uncertainty which was already taken into account in the low-fidelity step. The multifidelity solution is then computed re-centering the Taylor expansions in $M_{LF}(t_{k})$ on the high-fidelity trajectories $\{{\bf\it x}^{(i)}_{HF}(t_{k})\}$ . The $i^{th}$ domain of $M_{MF}(t_{k})$ is thus given by

[TABLE]

which is a Taylor polynomial centered at ${\bf\it x}^{(i)}_{HF}(t_{k})$ and whose non-constant coefficients model the uncertainty due to the ICs only. Information on the cumulative effects of stochastic accelerations is kept separately in $M_{PN}(t_{k})$ , and added to Eq. 41 when computing the polynomial bounds in observables space as described in Section 3.3. These two manifolds are collectively referred as $M_{{\bf\it x}}(t_{k})=\{([{\bf\it x}_{MF}^{(i)}(t_{k})],{\bf\it P}_{PN}^{(i)}(t_{k}))\}$ and constitute all the information on the SO state that is retained between subsequent steps.

For computational purposes, alternate equinoctial elements are used for propagation. This set of coordinates is defined in terms of the six Keplerian parameters as [32]

[TABLE]

with $a$ the semi-major axis, $e$ the eccentricity, $i$ the inclination, $\Omega$ the right ascension of the ascending node, $\omega$ the argument of periapsis, and $M$ the mean anomaly, where $n$ is the mean motion, and $\lambda$ the mean longitude. Given the quasi-linearity of the ODEs governing their evolution over time in perturbed Keplerian dynamics, this choice leads to smaller number of domains generated in the low-fidelity step compared to either Cartesian, Keplerian or (modified) equinoctial parameters [33, 34].

3.3 Projection onto observables space

Consider a set of osculating orbital elements ${\bf\it x}(t_{k})$ at epoch $t_{k}$ and an optical telescope as described in Section 2.1.2. The spherical coordinates $(\rho_{k},\alpha_{k},\delta_{k})$ of the SO as seen from the telescope at epoch $t_{k}$ are given by Eq. 6 with ${\bf\it\rho}=[\rho_{x}\ \rho_{y}\ \rho_{z}]^{T}$ LOS vector obtained from Eq. 4. Moreover, the SO’s position ${\bf\it r}_{k}$ is easily obtained by transforming ${\bf\it x}(t_{k})$ to Cartesian parameters. These operations are carried out in the DA framework using the LOADS algorithm to accurately map the uncertainties in ${\bf\it x}(t_{k})$ through the sequence of nonlinear transformations. However, the manifold $M_{MF}(t_{k})=\{[{\bf\it x}_{MF}^{(i)}(t_{k})]\}$ output of Section 3.2 cannot be directly processed if the effects of process noise on the projected states have to be taken into account. The last are in fact represented by the covariance matrices in $M_{PN}(t_{k})$ , and an intermediate manifold must be built to include this information. Note that each polynomial in $M_{MF}(t_{k})$ is function of the independent DA variables $(\delta{\bf\it\alpha},\delta{\bf\it\delta})$ while there is no concept of Taylor expansion in ${\bf\it P}_{PN}^{(i)}(t_{k})$ . Six additional independent variables $\delta{\bf\it x}_{k}$ are thus considered to describe ${\bf\it P}_{PN}^{(i)}(t_{k})$ in the DA framework as

[TABLE]

with $c$ the $z$ -score introduced in Section 2.1.2, ${\bf\it P}^{(i)}_{PN}(t_{k})={\bf\it V}^{(i)}{\bf\it\Lambda}^{(i)}\left({\bf\it V}^{(i)}\right)^{T}$ the covariance matrix eigendecomposition, and ${\bf\it\lambda}^{(i)}$ the main diagonal of ${\bf\it\Lambda}^{(i)}$ . For each element in $M_{{\bf\it x}}(t_{k})=\{([{\bf\it x}_{MF}^{(i)}(t_{k})],{\bf\it P}_{PN}^{(i)}(t_{k}))\}$ a Taylor expansion that represents the uncertainty due to both the ICs and the stochastic accelerations is thus given by

[TABLE]

which is a polynomial in twelve independent DA variables. These polynomials are then collected into a manifold $M_{MF,PN}(t_{k})$ which is subsequently projected onto the observables space $(\alpha,\delta)$ . The result is a manifold of expected measurements $M_{\hat{{\bf\it y}}}(t_{k})$ with entries

[TABLE]

note that the cardinality of $M_{\hat{{\bf\it y}}}(t_{k})$ is generally greater that that of $M_{MF,PN}(t_{k})$ since several splits might occur while evaluating the transformations in the DA framework. The domains $[{\bf\it x}_{MF,PN}^{(i)}(t_{k})]$ that maps onto $[\hat{{\bf\it y}}^{(i)}(t_{k})]$ are then easily reconstructed a posteriori to obtain a bijection between the elements of the two manifolds. This step is essential to correctly perform the subsequent domain pruning.

3.4 Domain pruning

Consider the manifold of observables $M_{\hat{{\bf\it y}}}(t_{k})$ estimated from $M_{MF,PN}(t_{k})$ as in Section 3.3 and the real measurements $(\alpha_{k},\delta_{k})$ taken from the same telescope at epoch $t_{k}$ . These measurements are assumed to be corrupted by uncorrelated white noise with standard deviations $\sigma_{\alpha_{k}},\sigma_{\delta_{k}}$ respectively. The pruning algorithm seeks to reduce the size of $M_{\hat{{\bf\it y}}}(t_{k})$ exploiting the new information obtained from the sensor. The real measurements are firstly initialized as a DA vector

[TABLE]

with $c$ the $z$ -score as in Eq. 7 and $\delta{\bf\it y}$ the first-order deviations in $\alpha_{k},\delta_{k}$ . Rigorous bounds on $[{\bf\it y}]$ are then given by the square interval $[\alpha_{k}\pm c\cdot\sigma_{\alpha_{k}}]\times[\delta_{k}\pm c\cdot\sigma_{\delta_{k}}]$ which is a simplification of Eq. 9 for first-order univariate polynomials.

For each domain in $M_{\hat{{\bf\it y}}}(t_{k})$ , its bounds are then estimated via Eq. 9 and compared with the ones above. If the two intervals overlap, the domain is retained for the subsequent steps. The polynomial is otherwise discarded from $M_{\hat{{\bf\it y}}}(t_{k})$ since it is unlikely to include the true, yet unknown, orbit state. Similarly, if no intersection is found between $M_{\hat{{\bf\it y}}}(t_{k})$ and $[{\bf\it y}]$ , the measurement is deemed to be an outlier and no pruning is performed. The observation is subsequently excluded from the set processed with the batch OD algorithm described in Section 2.3.1. As last step, pruning is also performed on $M_{MF}(t_{k})$ and $M_{PN}(t_{k})$ so as to reestablish a bijection between the elements of the three sets before processing the next measurement.

3.5 Combined prediction and pruning

The overall algorithm for the prediction and pruning of orbit states is presented in this section. Starting from the ICs given by Eq. 39, a sequential procedure is implemented to progressively tighten the bounds on $(\delta{\bf\it\alpha},\delta{\bf\it\delta})$ by exploiting the information available from subsequent measurements.

The estimated Cartesian state is firstly transformed into alternate equinoctial elements. This transformation is evaluated in the DA framework using the LOADS scheme and the transformed ICs are thus represented as a manifold of polynomials $M_{MF}(t_{0})$ function of $(\delta{\bf\it\alpha},\delta{\bf\it\delta})$ . The state covariances due to process noise are also initialized as $M_{PN}(t_{0})=\{{\bf\it 0}_{6\times 6}\}$ for each polynomial in $M_{MF}(t_{0})$ , and the two manifolds grouped into $M_{{\bf\it x}}(t_{0})=\{([{\bf\it x}^{(i)}(t_{0})],{\bf\it 0}_{6\times 6}^{(i)})\}$ .

The measurement collected at time $t_{k}$ is denoted by $(\alpha_{k},\delta_{k})$ , $1\leq k\leq M$ , and the associated $1\sigma$ uncertainties by $(\sigma_{\alpha_{k}},\sigma_{\delta_{k}})$ , such that the developed algorithm requires $M$ sequential steps to estimate $M_{{\bf\it x}}(t_{M})$ from $M_{{\bf\it x}}(t_{0})$ while exploiting the information available from the measurements. Each step $k$ is further subdivided into the following operations

Propagation of $M_{{\bf\it x}}(t_{k-1})$ to the next epoch $t_{k}$ using the multifidelity scheme in Section 3.2 to obtain $M_{{\bf\it x}}(t_{k})$ 2. 2.

Construction of the inflated manifold $M_{MF,PN}(t_{k})$ and projection of the last onto the observables space $(\rho,\alpha,\delta)$ as in Section 3.3 to obtain $M_{\hat{{\bf\it y}}}(t_{k})$ 3. 3.

Pruning of $M_{\hat{{\bf\it y}}}(t_{k})$ and $M_{{\bf\it x}}(t_{k})$ with information from $(\alpha_{k},\delta_{k})$ as in Section 3.4 to obtain the reduced manifolds $M^{\prime}_{\hat{{\bf\it y}}}(t_{k}),M^{\prime}_{{\bf\it x}}(t_{k})$ as well as isolate potential outliers. If an outlier is found, no pruning is performed and $M^{\prime}_{{\bf\it x}}(t_{k})$ coincides with $M_{{\bf\it x}}(t_{k})$ 4. 4.

Merging of the domains in $M^{\prime}_{{\bf\it x}}(t_{k})$ as in Section 2.2.4 since nonlinearities are weaker in orbital elements space than in $(\rho,\alpha,\delta)$ and fewer domains are thus required to accurately capture the uncertainty in this space

The procedure thus ends at epoch $t_{M}$ with the processing of the last available measurement. Two lists of indexes $K_{c}$ and $K_{o}$ corresponding to correlated and outlier measurements are also available at this stage.

3.6 Differential algebra orbit determination

Outputs of Section 3.5 include a manifold of osculating elements $M^{\prime}_{{\bf\it x}}(t_{M})$ at the last measurement epoch $t_{M}$ and a set of correlated measurements $({\bf\it y}_{k})_{k\in K_{c}}$ . From these data a batch OD problem can be then set up as follows.

Firstly, the domains $M^{\prime}_{{\bf\it x}}(t_{0})$ that map onto the pruned manifold at $t_{M}$ are reconstructed by applying the splitting history of $M^{\prime}_{{\bf\it x}}(t_{M})$ to the ICs $M_{{\bf\it x}}(t_{0})$ . Among all domains in $M^{\prime}_{{\bf\it x}}(t_{0})$ , the center of the one characterized by the smallest residuals with respect to the correlated observations $({\bf\it y}_{k})_{k\in K_{c}}$ is used as initial guess for the subsequent batch OD. The last is denoted as $\hat{{\bf\it x}}_{0}^{0}$ where the superscript $``0"$ refers to the zeroth iteration of the algorithm. The best estimate $\hat{{\bf\it x}}_{0}$ and associated covariance ${\bf\it P}_{0}$ are then computed using either the LS or LSAR algorithms described in Sections 2.3.2 and 2.3.3 respectively. Within each iteration $j$ of the OD scheme, DA is again leveraged to automate the computation of the observation function ${\bf\it h}_{k}$ defined in Eq. 21 and its partials ${\bf\it H}_{k}$ with respect to the current estimate $\hat{{\bf\it x}}_{0}^{j}$ . To do so, $\hat{{\bf\it x}}_{0}^{j}$ is firstly initialized as the DA vector $[\hat{{\bf\it x}}_{0}^{j}]=\hat{{\bf\it x}}_{0}^{j}+\delta{\bf\it x}_{0}$ . The propagation to $t_{k}$ and projection onto the observables space are then carried out in the DA framework after setting the expansion order equal to one. The Taylor expansion of the theoretical measurements $\hat{{\bf\it y}}_{k}$ at epoch $t_{k}$ is thus available as

[TABLE]

${\bf\it h}_{k}(\hat{{\bf\it x}}_{k})$ and ${\bf\it H}_{k}$ are then retrieved as the constant part and the linear part of Eq. 47, respectively. This avoids the needs for an analytical expression of $\partial{\bf\it h}_{k}/\partial{\bf\it x}_{k}$ and the integration of the $6\times 6$ variational equations for ${\bf\it\Phi}$ .

At convergence, the solution to the OD problem consists in an estimated state $\hat{{\bf\it x}}_{0}$ that minimizes the cost function Eq. 26 or Eq. 29 and an associated covariance matrix $\mathbf{P}_{0}$ . If a new set of observations become available, it is then conceivable to repeat the process and compute an updated estimate $\hat{{\bf\it x}}^{\prime}_{0}$ and covariance ${\bf\it P}^{\prime}_{0}$ by taking into account the newly available information. Rather than repeating each step in Fig. 1 by considering the two batches of measurements at once, it is more efficient to initialize $M_{{\bf\it x}}(t_{0})$ from $\hat{{\bf\it x}}_{0},{\bf\it P}_{0}$ and run the pruning and OD algorithms only on the second batch. Assuming a Gaussian distribution for $\hat{{\bf\it x}}_{0}$ , the last is initialized as the DA vector

[TABLE]

with $c$ the $z$ -score, ${\bf\it P}_{0}={\bf\it V}_{0}{\bf\it\Lambda}_{0}{\bf\it V}_{0}^{T}$ the covariance matrix eigendecomposition and ${\bf\it\lambda}_{0}$ the main diagonal of ${\bf\it\Lambda}_{0}$ . In this case $M_{{\bf\it x}}(t_{0})$ will contain a single element, i.e. $[\mathbf{x}_{0}]$ , which is a multivariate polynomial in the six independent DA variables $\delta{\bf\it x}_{0}$ . It is thus different from Section 3.5 where the initial manifold at $t_{0}$ is function of the angular variables $(\delta{\bf\it\alpha},\delta{\bf\it\delta})$ introduced in Section 3.1 and might contain more than one element if splits occurred during IOD. Moreover, representing the state estimate with a mean vector $\hat{{\bf\it x}}_{0}$ and a covariance matrix ${\bf\it P}_{0}$ is a limitation of the current algorithm, since all second order information gained with the DA-based IOD and pruning routines is lost when the batch OD algorithm is started. Future developments will thus focus on replacing the current OD routines with high-order DA-based sequential filters [35]. The pruning and OD steps could then be combined in a single sequential algorithm in which the state update is also performed in the algebra of Taylor polynomials. The final output would then be a second-order manifold $M_{{\bf\it x}}(t_{M})$ , and newly available observations for $t>t_{M}$ would be seamlessly processed without the needs for a re-initialization of the DA variables as in Eq. 48.

4 Numerical Application

This section presents an application of the developed algorithm in real operational scenarios. The selected object is the rocket body of an Ariane 44L launched in 2002 and catalogued with NORAD ID 27402. This object is on a geostationary transfer orbit (GTO) with semi-major axis of about $23\,000\text{\,}\mathrm{km}$ and inclination of $3.9\text{\,}\mathrm{\SIUnitSymbolDegree}$ , and is regularly tracked by the optical telescopes of the TAROT network [36]. This scenario was selected since () are notoriously challenging for both IOD and UP tasks. Specifically for the object of interest, observations are available only in the proximity of the apogee where the slow motion of the SO, coupled with a large distance from the sensors, leads to small angular separations between measurements and thus large uncertainties in the computed IOD solution. GTOs are then characterized by very fast perigee passages, during which the objects is subject to strong perturbations from the Earth’s atmosphere, and slow apogee passages, during which the effects of atmospheric drag are negligible, while Solar radiation pressure and lunisolar perturbations prevail. This continuous alternation of fast and slow dynamics, coupled with a large eccentricity that magnifies its nonlinearities, make more difficult to maintain an accurate representation of the state over time. A successful application of the algorithm to this scenario will thus give confidence on its robustness and applicability to a wide spectrum of orbital regimes.

Firstly, the performance of DA-based UP methods is assessed in Section 4.1 by comparing their accuracy and computational efficiency against a high-fidelity MC simulation. Two analyses are then conducted in Sections 4.2 and 4.3 to highlight different aspects of the proposed tool. In both cases, two objects are considered: the one of interest (the target), and a third-party one (the outlier). The outlier is simulated, and noised measurements are produced out of it to simulate observations that should not be associated to the target object. Instead, the observations that should be associated to the SO of interest are either generated from a synthetic trajectory in Section 4.2 or drawn from real observations in Section 4.3. Both analyses consider a time span of about five days that includes five passages of the SOs over the two telescopes in La Reunion and Calern (France). Passages epochs, duration and number of measurements are summarized in Table 1.

In the followings the analytical SGP4 propagator [30] is used as low-fidelity dynamical model. Instead, a high-fidelity numerical propagator is set up with the force models listed below

•

Earth’s gravity potential modeled with spherical harmonics truncated at order/degree 16 [37]

•

Sun and Moon third-body attraction with bodies’ position from JPL DE440 ephemerides [38]

•

Atmospheric drag with Harris-Priester atmosphere model [39]

•

Solar radiation pressure with cylindrical Earth’s shadow model [40]

Taylor polynomials in Eqs. 35, 43 and 46 are initialized with a $z$ -score equal to $c=3$ such that the $99.73\%$ of the probability mass is accurately mapped through each nonlinear transformation using the LOADS scheme. The algorithms were coded in Java and make use of the following external libraries

•

CNES’s PACE (Polynomial Algebra Computational Engine) library as DA backend to perform operations on Taylor polynomials

•

Orekit333https://www.orekit.org/ for standard flight dynamics routines (SGP4 propagator, force models, time scales and reference frames conversions)

•

Hipparchus444https://hipparchus.org/ for the 8(5,3) Dormand-Prince numerical integrator and the Levenberg–Marquardt optimizer

•

MOSEK555https://www.mosek.com/ interfaced through its Fusion API for Java666https://docs.mosek.com/latest/javafusion/index.html as LPP solver

All simulations were run on a cluster equipped with Skylake Intel® Xeon® Gold 6126 @ $2.6\text{\,}\mathrm{GHz}$ and $96\text{\,}\mathrm{GB}$ of RAM running Red Hat® Enterprise Linux® 7.8.

4.1 Validation of multifidelity method

The performance of the multifidelity UP method presented in Section 3.2 is discussed in this section to justify its adoption within the developed OD framework. Four different UP schemes are considered, namely

Low-fidelity (LF) DA-based propagation (first step of Section 3.2) 2. 2.

Multifidelity (MF) DA-based propagation (as used within the OD framework) 3. 3.

High-fidelity (HF) DA-based propagation (DA-aware numerical propagator directly embedded in the LOADS framework) 4. 4.

MC simulation with $N_{s}=10^{4}$ samples using the high-fidelity dynamics

A comparison is then performed in terms of both accuracy and computational efficiency, with the MC simulation taken as ground truth. A manifold of ICs $M(t_{0})=\{[{\bf\it x}^{(i)}(t_{0})]\}$ is firstly obtained by solving the IOD problem as described in Section 3.1 using the real measurements from the first passage in Table 1. These ICs are then propagated over a five days time span with the four methods listed above. First, random samples $\{{\bf\it x}^{(l)}(t_{0})\}$ for $l=1,\ldots,N_{s}$ are drawn from $M(t_{0})$ and propagated point-wise to obtain the reference empirical distribution $\{{\bf\it x}^{(l)}(t_{f})\}$ at final time. Then, $M(t_{0})$ is propagated in the DA/LOADS framework to obtain the three manifolds $M_{LF}(t_{f}),M_{MF}(t_{f})$ and $M_{HF}(t_{f})$ corresponding to the three DA-based UP schemes. These manifolds are then evaluated on $\{{\bf\it x}^{(l)}(t_{0})\}$ to obtain $\{{\bf\it x}_{LF}^{(l)}(t_{f})\},\{{\bf\it x}_{MF}^{(l)}(t_{f})\}$ and $\{{\bf\it x}_{HF}^{(l)}(t_{f})\}$ without the needs for further propagation. The root mean square error (RMSE) is finally used to assess the accuracy of the polynomial-based UP methods. The last is defined as

[TABLE]

with $\{{\bf\it x}^{(l)}\}$ the expected samples, i.e. the output of the MC simulation, and $\{\hat{{\bf\it x}}^{(l)}\}$ the actual samples, i.e. the result of the polynomial evaluation of either of the three manifolds at final time.

The RMSEs for the case under consideration are reported in Table 2. As seen from this table, the low-fidelity dynamics is unable to accurately capture the time evolution of the state distribution, resulting in average errors that exceed $218\text{\,}\mathrm{km}$ in position and $3.36\text{\times}{10}^{-2}\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ in velocity. The results are greatly improved with the introduction of the multifidelity correction, with RMSEs of about $2.24\text{\,}\mathrm{km}$ and $3.18\text{\times}{10}^{-4}\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ respectively. The last are further reduced with the high-fidelity DA-based propagation whose accuracy is in the order of $1.74\text{\times}{10}^{-1}\text{\,}\mathrm{km}$ in position and $2.41\text{\times}{10}^{-5}\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ in velocity, an order of magnitude better than the multifidelity solution.

However, as shown in Table 3, additional accuracy comes at a high computational cost. The high-fidelity propagator is in fact 43 times slower than its multifidelity counterpart, which is in turn 14 times slower than the low-fidelity one. Recalling that the last coincides with the first step of the multifidelity scheme, it is possible to conclude that the computational cost of the hybrid method is driven by the correction step, i.e. the point-wise propagation of the polynomials’ constant parts in high-fidelity dynamics. The overall cost of the multifidelity scheme is thus comparable to that of a high-fidelity MC simulation with a number of samples equal to the number of domains generated by the LOADS algorithm. In this case no splits are triggered during the propagation in alternate equinoctial elements, and a single split is required to map the uncertainty onto Cartesian space. The runtime of the hybrid scheme is thus similar to that of the LF simulation plus the time required to propagate a single MC sample in high-fidelity. This fact can be easily verified by comparing the first, second and fourth columns of Table 3.

Figures 2 and 3 depicts the MC samples (gray plus signs) and the edges of the polynomial domains for the four cases enumerated before. Domains in blue correspond to the low-fidelity solution, which is clearly shifted with respect to the ground truth in both $y-z$ and $v_{y}-v_{z}$ projections. Instead, domains in gold (multifidelity solution) and in black (high-fidelity solution) overlap almost perfectly and provide a visual demonstration of the accuracy of the proposed multifidelity propagation algorithm. As expected, most of the MC samples fall inside these edges thus demonstrating the validity of the DA-based UP methods. Moreover, the good accuracy of the multifidelity method can be explained by examining the different steps involved to obtain the final solution. The first one consists in propagating the initial uncertainty in low-fidelity, and the output is equivalent to the low-fidelity solution depicted in blue. As seen in both pictures, the union of the regions mapped by these polynomials has the same shape of that mapped by the high-fidelity ones (delimited by black edges), and the only noticeable difference is a shift between the two. It is thus possible to conclude that SGP4, even though it fails in capturing the evolution of the absolute dynamics, is capable of accurately capture the relative motion in the neighborhoods of the polynomials’ constant parts delimited by the depicted edges. The second step is instead a propagation of the polynomials’ constant part in high-fidelity using the same dynamical model as the MC simulation. Once these two solutions are available, the high-fidelity trajectories are taken as reference with respect to which the uncertainty is expressed, and a simple substitution of the polynomials’ centers allows to obtain a solution that matches the accuracy of the high-fidelity one. To conclude, the multifidelity solution provides the best trade-off for the use case under investigation. The small loss in accuracy observed in Table 2 is in fact compensated by the large computational gain shown in Table 3.

4.2 Synthetic measurements

The first analysis of the complete OD chain uses synthetic measurements for both objects, so as a true reference trajectory is available to assess the algorithm’s performance. For the target SO an assumed true state is retrieved from TLE data. The outlier orbit is then obtained by perturbing the target’s eccentricity by $\Delta e=-0.02$ . The Keplerian parameters at epoch are given in Table 4 for both SOs. The trajectories of the two objects are simulated using the high-fidelity propagator and synthetic measurements are drawn from them. These are then corrupted with a zero-mean Gaussian white noise with standard deviations $\sigma_{\alpha}=$ 1.285\text{,}\mathrm{\SIUnitSymbolArcsecond} $and $\sigma_{\delta}=$1.280\text{\,}\mathrm{\SIUnitSymbolArcsecond}$ to simulate sensor noise. Four different scenarios are then analyzed as listed in Table 5. In scenario A and B only correlated measurements are processed. These are meant to validate the tools before further complexity is added to the problem. Observations associated to the outlier are then added in scenarios C and D.

The IOD algorithm described in Section 3.1 is run on the first passage and the estimated state at $t_{0}$ used to initialize the pruning scheme. In this case, the LOADS algorithm generates a single domain and thus $M_{{\bf\it x}}=\{[{\bf\it x}_{0}]\}$ . In scenarios B and D this manifold is then propagated to the latest epoch using the algorithm in Section 3.5. The evolution over time of the number of domains is summarized in Tables 6 and 7 respectively. A visual representation is also given in Figs. 4 and 5. In this pictures, orange squares and yellow crosses correspond to the size of the osculating elements manifolds after propagation (step #1) and merging (step #4). Blue diamonds and green plus signs correspond to the size of the observables manifolds after projection (step #2) and pruning (step #3). Vertical dashed lines correspond to the observations epochs.

The effectiveness of the pruning algorithm is clearly seen by comparing, for the same epoch, the number of domains after steps #2 and #3. It is also noted that the merging step is key in limiting the number of domains processed in step #1. If not performed, the propagation algorithm will enter the next iteration with a number of domains corresponding to the third column (or green plus signs) rather than the fourth one (or yellow crosses). The effectiveness of the selected state representation in limiting the number of domains generated during propagation is also demonstrated. All splits are in fact triggered by the projection onto the observables space and the increased number of domains to be propagated in subsequent steps is the result of an imperfect merging during step #4. Since domains are discarded in step #3, several incomplete triplets are found by the merging scheme which thus fails to recombine the polynomials as they were before step #2. The correct detection of outliers is then demonstrated in Table 7, where no modification of the domains occur in the third measurement block which corresponds to uncorrelated measurements. In this case, the algorithm correctly detects the outliers and does not prune the projected manifold.

The success rate of the pruning algorithm in correctly identifying outlier measurements is reported in Table 8 for scenarios without outliers (cases A and B) and with outliers (cases C and D) respectively. This table shows the number of correctly associated (true positives) and discarded (true negative) measurements as well as the number of those that should have been marked as correlated (false negatives) or outliers (false positives) but they have been not. Results obtained with both low-fidelity and multifidelity UP methods are shown. The high-fidelity solution is instead omitted since a perfect match is already obtained with the multifidelity one. The same is however not true for the low-fidelity method for which five false negatives are observed in both scenarios. The importance of the multifidelity correction is thus once again remarked by these data.

An example of pruning is presented in Fig. 6 where the estimated bounds for both $M_{\hat{{\bf\it y}}}(t_{k})$ (gray boxes) and $M_{\hat{{\bf\it y}}}^{\prime}(t_{k})$ (gold boxes) are plotted in $(\alpha,\delta)$ space. The large overlap between gold boxes is an indicator that most of the uncertainty is on the LOS distance which cannot be observed with optical telescopes. Still, the pruning algorithm is highly effective in reducing the uncertainty on the observed quantities. This is demonstrated also in Fig. 7 which represents the uncertainty on the IOD solution in $(\delta_{2},\delta_{3})$ space. The gray box identifies the edges of the initial domain which corresponds to the uncertainty on the measurements. Instead, the golden boxes are the edges of the domains that maps into the solution manifold at the last epoch. It is thus seen the effectiveness of the algorithm in tightening the bounds on the initial estimate as the propagation steps forward and more information is exploited to refine the IOD solution.

OD is then performed for all scenarios comparing the performance of the LS and LSAR algorithms introduced in Section 2.3.1. The true Cartesian state of the target object and the IOD solution used to initialize the OD routines are given in Table 9. The resulting estimation errors and $\pm 3\sigma$ bounds are summarized in Tables 10, 11 and 12 for scenarios A and B, C and D respectively.

Table 10 corresponds to the two scenarios with no outliers in which the pruning algorithm has no effect on the final solution. Both algorithms provide a consistent estimate of the SO state converging in three and two iterations for the LS and LSAR methods respectively. As expected, the LS solution is characterized by tighter bounds on the states since the algorithm is the minimum variance estimator for a model with zero-mean Gaussian noise distribution [15].

The effects of measurement outliers on the estimated state is then shown in Table 11. With no measurement preprocessing, the LS estimator fails to converge to the true solution and returns an estimate characterized by large errors in both position ( $\approx$ 700\text{,}\mathrm{km} $) and velocity ($\approx$0.12\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ ) converging in six iterations. At the same time, the variance on the estimated state is close to that of scenarios A and B, thus placing the true state outside the $\pm 3\sigma$ bounds. On the contrary, the LSAR algorithm is capable of rejecting the outliers and still converging to the true solution at the expenses of a longer runtime (five iterations vs two in the previous cases). The estimated bounds are also consistent to those obtained in scenarios A and B and they thus include the true state reported in Table 9.

The accuracy of the LS estimate is then restored in scenario D in which outliers have been detected and discarded by the pruning scheme. Results in Table 12 are in fact similar to those in Table 10 and the LS estimation is again characterized by tighter bounds on the states. The LSAR algorithm also benefits from measurements preprocessing since excluding the uncorrelated observations from the processed ones reduces the number of iterations from five to just two. This results in a total runtime of about 40% of that of scenario C as can be verified by inspecting in Table 13.

4.3 TAROT measurements

The second analysis uses real observation data from the TAROT network to estimate the target SO state. Synthetic measurements for the outlier orbit generated in Section 4.2 are still used to corrupt the available data and assess the robustness of the proposed tool. The four scenarios presented in Section 4.2 are reproduced here as summarized in Table 14. To account for the unknown forces acting on the SOs, stochastic accelerations modeled as additive white Gaussian noise are included in the high-fidelity dynamics. Their effect is taken into account by the pruning scheme using the SNC algorithm described in Section 2.1.1. A deterministic dynamics is instead used for the subsequent batch OD.

The IOD algorithm run on the first passage still return a single domain which is then used to initialize the pruning scheme in scenarios F and H. The time evolution of the number of domains is summarized in Tables 15 and 16 and depicted in Figs. 8 and 9 respectively for the two scenarios. The immediate effect of stochastic accelerations can be identified in the larger number of domains generated during the projection (#2) step. Indeed, integrating the uncertainty due to these perturbations results in an inflated covariance for the prior states which, in turns, triggers more splits when projected onto the observables space. At the same time, the pruning scheme is even more effective and the ratio between projected and retained domains is higher in Table 15 than in Table 6. The success rate of the pruning algorithm in correctly identify outlier measurements is summarized in Table 17 and is identical to that of Table 8. The effectiveness of the multifidelity UP method in providing an accurate estimate of the prior at a fraction of the computational cost of a high-fidelity propagation is thus demonstrated also in this case.

Two examples of pruning are depicted in Figs. 10 and 11 for scenarios F and H respectively. These correspond to the epochs with the highest ratio of projected over retained domains, i.e. when the largest number of domains is dropped and the uncertainty is reduced the most. As expected, for scenario H this occurs at the first measurements after the block of outliers which ends a long propagation arc without useful information to reduce the growing uncertainty (since uncorrelated measurements are ignored).

The uncertainty on the IOD solution for the two scenarios is depicted in Figs. 12 and 13 in which a projection of the initial (gray box) and retained domains (golden boxes) onto the $(\delta_{2},\delta_{3})$ space is shown. A comparison of the area covered by the golden boxes demonstrates that for both scenarios the correct subregion of the initial domain that contains the true measurement is identified. The last corresponds to the lower-right portion of the initial rectangular domain, i.e. to above-average values of $\delta_{2}$ and below-average values of $\delta_{3}$ .

Similarly to Section 4.2, OD is performed in the four scenarios to assess the performance of the LS and LSAR estimators while processing both raw and retained observations. The true state is not available in this case, and the best estimate available from the CNES’s catalog is used as ground truth. The last was obtained by processing measurements spanning several years of observations, and it is continuously updated as soon as new data is available using a very high-fidelity dynamical model to propagate the best estimate to each measurement epoch. For the shake of limiting the overall runtime, this paper considers only a very limited number of observations, and no prior information is assumed except that all measurements are correlated to the same SO. The high-fidelity dynamical model was also simplified for the same reason. As a consequence, errors in the order of $10\text{\,}\mathrm{km}$ and $0.001\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ over a five days time span are considered acceptable in this case. The aforementioned reference solution is thus reported in Table 18 together with the initial guesses provided to the batch estimators. The lasts coincide with the raw IOD solution for scenarios E and G (since no pruning is performed) and to the domain that is characterized by the smallest residuals among the retained ones for scenarios F and H.

Table 19 reports the LS and LSAR estimation errors and $\pm 3\sigma$ bounds for scenario E. These are qualitatively identical to those of scenario F which are omitted for conciseness. The LS solution is again characterized by tighter bounds on the states, even though the errors with respect to the assumed reference are larger than for the LSAR estimate. In this case, three iterations are needed for both estimators to converge.

The effects of measurements outliers among the processed set is then shown in Table 20 which summarizes the results obtained for scenario G. In this case, five and three iterations are needed by the LS and LSAR algorithms to converge. As for scenario C, the LS estimate exhibits large errors in both positions ( $\approx$ 550\text{,}\mathrm{km} $) and velocity ($\approx$0.1\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ ) thus supporting the needs for a preprocessing algorithm capable of detecting and isolating uncorrelated observations. On the contrary, the LSAR scheme shows an excellent performance that matches the one of scenario H reported in Table 21. At the same time, the presence of uncorrelated measurements among the processed set does not influence the estimated variance as can be verified by comparing the second and fourth rows in Tables 19 and 20 respectively.

Looking at Table 22, no reduction in the number of iterations is observed when processing correlated measurements only instead of the entire available set. The inclusion of the pruning step is however paramount to restore the accuracy of the LS estimate as demonstrated in Table 21. In this case a consistent estimate is in fact obtained and the LS estimator converges to a solution closer to the true one than that of Table 19.

5 Conclusions

A pruning algorithm for the refinement of an IOD solution and the detection of measurement outliers was presented. This scheme is broken down into four main steps, carried out sequentially for each available observation: propagation of the state estimate, projection onto the observables space, domains pruning based on the matching between predicted and actual measurements, and domains merging to minimize the number of trajectories to be further propagated. The uncertainty on the IOD solution is represented as a Taylor polynomial and DA techniques are leveraged to efficiently propagate the uncertainty over time. A novel splitting algorithm is also employed for the automatic control of the error committed by the truncated polynomials. To further reduce its computational effort, the prediction step is based on a bifidelity dynamical model that couples a DA-aware SGP4 propagator with an high-fidelity numerical integrator. If the first is used to map the initial uncertainty on the state to the next observation epoch, the last is employed to compute an accurate reference trajectory on which said uncertainty is centered and to include the effects of stochastic accelerations in the predicted estimate. After projecting the prior distribution onto the observables space, polynomial bounding techniques are used to estimate the range of variation of each domain and prune the ones that do not overlap with observed data at the same epoch, thus reducing the volume of the propagated uncertainty and tightening the bounds on the initial IOD solution. Measurements outliers are also detected in this step if no intersection between predicted and actual measurements is found. This algorithm is then employed as a preprocessing scheme to filter the information used by subsequent batch OD tools. Its performance is demonstrated for a GTO object using synthetic as well as real observation data from the TAROT network. It is verified to greatly improve the accuracy of LS estimates, notoriously affected by the presence of outliers among the processed observations, and to reduce the runtime of LSAR estimators, which converge in less iterations if processing correlated measurements only instead of the entire available set. The successful application of this method to real data in a challenging scenario, namely a object characterized by poor observability and subject to a strongly nonlinear dynamics, suggests that it is applicable to the full spectrum of orbit regimes.

Acknowledgments

This work is co-funded by the CNES through A. Fossà PhD program, and made use of the CNES orbital propagation tools, including the PACE library. The authors gratefully acknowledge Dr. Laura Pirovano for her support during the implementation and validation of the IOD algorithm.

Bibliography40

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] Ryan S. Park and Daniel J. Scheeres “Nonlinear Mapping of Gaussian Statistics: Theory and Applications to Spacecraft Trajectory Design” In Journal of Guidance, Control, and Dynamics 29.6 American Institute of Aeronautics Astronautics Inc., 2006, pp. 1367–1375 DOI: 10.2514/1.20177 · doi ↗
2[2] Martin Berz “Modern Map Methods in Particle Beam Physics” London: Academic Press, 1999, pp. 317 URL: http://bt.pa.msu.edu/cgi-bin/display.pl?name=AIEP 108book
3[3] M. Valli, R. Armellin, P. Di Lizia and M.. Lavagna “Nonlinear Mapping of Uncertainties in Celestial Mechanics” In Journal of Guidance, Control, and Dynamics 36.1 , 2013, pp. 48–63 DOI: 10.2514/1.58068 · doi ↗
4[4] S.J. Julier and J.K. Uhlmann “Unscented Filtering and Nonlinear Estimation” In Proceedings of the IEEE 92.3 , 2004, pp. 401–422 DOI: 10.1109/JPROC.2003.823141 · doi ↗
5[5] Nagavenkat Adurthi and Puneet Singla “Conjugate Unscented Transformation-Based Approach for Accurate Conjunction Analysis” In Journal of Guidance, Control, and Dynamics 38.9 American Institute of Aeronautics Astronautics, 2015, pp. 1642–1658 DOI: 10.2514/1.G 001027 · doi ↗
6[6] Brandon A. Jones, Alireza Doostan and George H. Born “Nonlinear Propagation of Orbit Uncertainty Using Non-Intrusive Polynomial Chaos” In Journal of Guidance, Control, and Dynamics 36.2 American Institute of Aeronautics Astronautics Inc., 2013, pp. 430–444 DOI: 10.2514/1.57599 · doi ↗
7[7] Ya-zhong Luo and Zhen Yang “A review of uncertainty propagation in orbital mechanics” In Progress in Aerospace Sciences 89 Elsevier, 2017, pp. 23–39 DOI: 10.1016/j.paerosci.2016.12.002 · doi ↗
8[8] Alexander Wittig et al. “Propagation of large uncertainty sets in orbital dynamics by automatic domain splitting” In Celestial Mechanics and Dynamical Astronomy 122.3 Kluwer Academic Publishers, 2015, pp. 239–261 DOI: 10.1007/s 10569-015-9618-3 · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

A Multifidelity Approach to Robust Orbit Determination

Abstract

1 Introduction

2 Mathematical background

2.1 Modeling of the problem

2.1.1 State noise compensation

2.1.2 Measurement process and sensor noise

2.2 Differential algebra

2.2.1 Polynomial bounding

2.2.2 Uncertainty propagation

2.2.3 Low-order automatic domain splitting

2.2.4 Domain merging

2.3 Orbit determination

2.3.1 Batch orbit determination algorithms

2.3.2 Least squares

2.3.3 Minimum L1L_{1}L1​ norm estimates

3 Robust orbit determination

3.1 Initial orbit determination

3.2 Multifidelity uncertainty propagation

3.3 Projection onto observables space

3.4 Domain pruning

3.5 Combined prediction and pruning

3.6 Differential algebra orbit determination

4 Numerical Application

4.1 Validation of multifidelity method

4.2 Synthetic measurements

4.3 TAROT measurements

5 Conclusions

Acknowledgments

2.3.3 Minimum $L_{1}$ norm estimates