One possible approach is to try RANSAC. RANSAC is applicable where we expect: (a) some fraction of the observations (the rays) are erroneous (outliers) and should be ignored, and (b) the remaining ones are all approximately consistent with some low-dimensional probabilistic model. The RANSAC method then tries to find the best model that (a) minimizes the number of outliers, and (b) among all such models, is as consistent as possible with the remaining observations.
In your case, RANSAC might work well enough. Some small fraction of the rays are in a completely wrong direction, and will be ignored as outliers. The remaining ones should all be roughly consistent. The model could be (for instance) a two-dimensional Gaussian with diagonal covariance matrix. Such a model has 4 free parameters, so given 2 or more rays, you can fit a Gaussian to it. Given a model, you can identify which rays are inliers (plausibly consistent with the model) and which rays are outliers.
Given this, one way to apply RANSAC to your problem would be to repeat the following recipe 1000 times:
- randomly pick 5 rays
- fit a Gaussian to those 5 rays
- for each ray, compute the likelihood that it came from that Gaussian, and throw away all those whose likelihood is below some threshold (these are the outliers)
- fit a new Gaussian to all those remaining rays (the inliers)
- count how many rays remain consistent with the revised model (i.e., count the number of inliers)
Then keep the model with the largest number of inliers, and use the corresponding Gaussian as your model. In particular, the mean of that Gaussian is your best guess for the target, and the equal-likelihood contour of that Gaussian (the counter corresponding to 2 standard deviations away from the mean) is your circle. All of the above constants (1000, 5, 2) can be tuned, and variants of this algorithm can be considered.
Why does this work? Suppose 80% of the rays are correct and 20% are thrown in the completely wrong direction. Then a random sample of 5 of the rays will, with probability $0.8^5 \approx 0.33$, contain only good rays, and when you fit a Gaussian to them, will give a reasonable estimate of the target and the surrounding circle. The remaining 67% of samples will have at least one bad ray, which causes a bad model, but the bad models will be thrown away.
The general idea of this sort of approach is: (1) identify a scoring function, which given a model gives you a score for the model; (2) repeatedly pick a random sample of a few rays, fit a model to that sample, use the model to identify outliers, discard the outliers and fit a model to the remaining inliers, and then compute the score of the resulting candidate model; keep the best model you've ever seen. You can play with the type of model you allow, and the scoring function you use.
Some form of expectation-maximization (EM) potentially might yield even better results.