using System; using System.Collections.Generic; using System.Linq; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Security.Policy; using System.Text; using System.Threading.Tasks; using DW5S.Entity; namespace DW5S.KxcApi { /// ///TODO 定位算法接口.在这里调用罗博士的算法库 /// public static class PosApi { //置信度统计次数 private static int confidenceCount = 1000; //置信度有效距离(m) private static int confidenceDistance = 10000; #region cpp dll Interop //两星一地和三星解析定位算法(精度差、速度快,适用于置信度等统计时的多次调用) private const string gzPos = @"AddIns\定位\DLL_GZDW.dll"; //一星一地测向带参定位 private const string XDCX = @"AddIns\定位\DLL_DTO_DOA_DW.dll"; //三星双时差带参、三星双时差无参、三星双频差带参、双星时频差带参、两星一地无参定位及时差线、一星两地 private const string OtherPos = @"AddIns\定位\Position-New.dll";//DLL_11J_DW [DllImport(gzPos, EntryPoint = "DW_Analysis", CallingConvention = CallingConvention.Cdecl)]//两星一地和三星的解析定位算法(精度差、速度快) private extern static void DW_Analysis(double[] mainSatXYZ, double[] adja1XYZ, double[] adja2XYZ, double[] refStation, double tarDto1, double tarDto2, double refDto1, double refDto2, double[] posRes, int flag); [DllImport(XDCX, EntryPoint = "XD_CX_DW", CallingConvention = CallingConvention.Cdecl)]//一星一地测向带参 private extern static void X1D1_Pos20240305_Core(double[] mainSat, double[] satStation, double[] cdbStation, double[] cxStation, double[] refStation, double[] zone, double theta, double tarDto, double refDto, double[] res); [DllImport(OtherPos, EntryPoint = "SC_2X1D_DW", CallingConvention = CallingConvention.Cdecl)]//两星一地带参 private extern static void X2D1_Pos20240305_Core(double[] mainSat, double[] adjaSat, double[] cdbStation, double[] satStation1, double[] satStation2, double[] satStation3, double[] satStation4, double[] satStation5, double[] refStation, double[] zone, double tarSxDto, double tarXdDto, double samp_main_dto, double samp_neigh_dto, double[] res); [DllImport(OtherPos, EntryPoint = "XingDi_2X1D_DW_NoRef", CallingConvention = CallingConvention.Cdecl)]//两星一地无参 private extern static void X2D1_PosNoRef20240305_Core(double[] mainSat, double[] adjaSat, double[] cdbStation, double[] mainSatTarStation, double[] adjaSat1TarStation, double[] adjaSat2TarStation , double[] zone, double tarSxDto, double tarXdDto, double[] res); [DllImport(OtherPos, EntryPoint = "SanXing_DW", CallingConvention = CallingConvention.Cdecl)]//三星双时差带参 private extern static void X3_Pos20240305_Core(double[] mainSat, double[] adjaSat1, double[] adjaSat2, double[] mainSatTarStation, double[] adjaSat1TarStation, double[] adjaSat2TarStation , double[] mainSatRefStation, double[] adjaSat1RefStation, double[] adjaSat2RefStation, double[] refStation, double[] zone, double tarDto1, double tarDto2, double refDto1, double refDto2, double[] res); [DllImport(OtherPos, EntryPoint = "SanXing_DW_NoRef", CallingConvention = CallingConvention.Cdecl)]//三星双时差无参 private extern static void X3_PosNoRef20240305_Core(double[] mainSat, double[] adjaSat1, double[] adjaSat2, double[] mainSatTarStation, double[] adjaSat1TarStation, double[] adjaSat2TarStation , double[] zone, double tarDto1, double tarDto2, double[] res); [DllImport(OtherPos, EntryPoint = "TriStar_2DFO_DW", CallingConvention = CallingConvention.Cdecl)]//三星双频差带参 private extern static void X3_PosTwoDfo20240305_Core(double[] mainSat, double[] adjaSat1, double[] adjaSat2, double[] mainSatTarStation, double[] adjaSat1TarStation, double[] adjaSat2TarStation , double[] refStation, double[] zone, double tarDfo1, double tarDfo2, double refDfo1, double refDfo2, double fu1, double fd1, double fu2, double fd2, double[] res); [DllImport(OtherPos, EntryPoint = "TwoStar_DTFO_DW", CallingConvention = CallingConvention.Cdecl)]//双星时频差带参 private extern static void X2_Pos20240305_Core(double[] mainSat, double[] adjaSat, double[] mainSatTarStation, double[] adjaSatTarStation, double[] refStation, double[] zone , double tarDto, double tarDfo, double refDto, double refDfo, double fu1, double fd1, double fu2, double fd2, double[] res); #endregion /// /// 一星一地带参,返回经度、纬度、高度、镜像点、置信度,数组长度为7 /// /// 参估结果 /// 站点信息 /// 测向结果 /// 是否计算置信度 /// public static double[] X1D1_Pos(CgRes cgRes, StationRes sRes, CxRes cxRes, bool CalcConfidence = false) { if (cgRes.Snr1 == 0) return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; if (cgRes.RefYbDto1 == null) return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; double[] mainSat = new double[3] { cgRes.MainX, cgRes.MainY, cgRes.MainZ }; double[] satStation = new double[3] { sRes.SatTxLon, sRes.SatTxLat, 0 }; double[] cdbStation = new double[3] { sRes.CdbTxLon.Value, sRes.CdbTxLat.Value, 0 }; double[] cxStation = new double[3] { sRes.CxLon.Value, sRes.CxLat.Value, 0 }; double[] refStation = new double[3] { sRes.RefLon.Value, sRes.RefLat.Value, 0 }; double dtoCdb = cgRes.Dto1 / 1e6; double[] zone = new double[] { -85, 85, -180, 180 }; //定位区域 double theta = cxRes.Fx;//单位° double[] res = new double[6]; var ybDto1 = cgRes.RefYbDto1.Value / 1e6; X1D1_Pos20240305_Core(mainSat, satStation, cdbStation, cxStation, refStation, zone, theta, dtoCdb, ybDto1, res); var posRes = ConvertToGeoPoint(res); if (CalcConfidence && IsGeoPoint2(posRes)) { int succeedCount = 0; for (int i = 0; i < confidenceCount; i++) { var mainSatNew = new double[3] { cgRes.MainX+RandomHelper.Normal(0,200), cgRes.MainY+RandomHelper.Normal(0,200), cgRes.MainZ+RandomHelper.Normal(0,200) }; var dtoCdbNew = (dtoCdb * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; var ybDtoNew = (ybDto1 * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; var thetaNew = theta + RandomHelper.Normal(0, 0.4);//单位° X1D1_Pos20240305_Core(mainSatNew, satStation, cdbStation, cxStation, refStation, zone, thetaNew, dtoCdbNew, ybDtoNew, res); var posResNew = ConvertToGeoPoint(res); if (IsGeoPoint2(posResNew)) { var distance = PhysicsHelper.DistanceGeo((posRes[0], posRes[1], 0), (posResNew[0], posResNew[1], 0)); if (distance <= confidenceDistance) { succeedCount++; } } } posRes[6] = (int)(succeedCount / (double)confidenceCount * 100); } return posRes; } /// /// 两星一地带参解析定位(精度差速度快,主要用于置信度分析等需要多次调用的情况) /// /// /// /// public static double[] X2D1_GzPos(CgRes cgRes, StationRes sRes) { if (cgRes.Snr1 == 0 || cgRes.Snr2 == 0) { return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; } double[] mainSat = new double[3] { cgRes.MainX, cgRes.MainY, cgRes.MainZ }; double[] adjaSat = new double[3] { cgRes.Adja1X.Value, cgRes.Adja1Y.Value, cgRes.Adja1Z.Value }; double[] cdbStation = new double[3] { sRes.CdbTxLon.Value, sRes.CdbTxLat.Value, 0 }; double[] refStation = new double[3] { sRes.RefLon.Value, sRes.RefLat.Value, 0 }; double dto1 = cgRes.Dto1 / 1e6; double dtoCdb = cgRes.Dto2.Value / 1e6; double ybDto1 = cgRes.RefYbDto1.Value / 1e6; double ybDto2 = cgRes.RefYbDto2.Value / 1e6; double ybDto = ybDto1 - ybDto2; double[] posRes = new double[3]; DW_Analysis(mainSat, adjaSat, cdbStation, refStation, dto1, -dtoCdb, ybDto, -ybDto1, posRes, 1); return posRes; } /// /// 两星一地带参,返回经度、纬度、高度、镜像点、置信度,数组长度为7 /// /// 参估结果 /// 站点信息 /// 是否计算置信度 /// public static double[] X2D1_Pos(CgRes cgRes, StationRes sRes, bool CalcConfidence = false) { if (cgRes.Snr1 == 0 || cgRes.Snr2 == 0) { return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; } if (cgRes.RefYbDto1 == null || cgRes.RefYbDto2 == null) { return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; } double[] mainSat = new double[3] { cgRes.MainX, cgRes.MainY, cgRes.MainZ }; double[] adjaSat = new double[3] { cgRes.Adja1X.Value, cgRes.Adja1Y.Value, cgRes.Adja1Z.Value }; double[] satStation = new double[3] { sRes.SatTxLon, sRes.SatTxLat, 0 }; double[] cdbStation = new double[3] { sRes.CdbTxLon.Value, sRes.CdbTxLat.Value, 0 }; double[] refStation = new double[3] { sRes.RefLon.Value, sRes.RefLat.Value, 0 }; double dto1 = cgRes.Dto1 / 1e6; double dtoCdb = cgRes.Dto2.Value / 1e6; double ybDto1 = cgRes.RefYbDto1.Value / 1e6; double ybDto2 = cgRes.RefYbDto2.Value / 1e6; double[] zone = new double[] { -85, 85, -180, 180 }; //定位区域 double[] res = new double[6]; X2D1_Pos20240305_Core(mainSat, adjaSat, cdbStation, satStation, satStation, satStation, satStation, satStation, refStation, zone, dto1, dtoCdb, ybDto1, ybDto2, res); ConvertToGeoPoint(res); var posRes = ConvertToGeoPoint(res);//精确搜索值 if (CalcConfidence && IsGeoPoint2(posRes)) { var posResGz = X2D1_GzPos(cgRes, sRes); if (IsGeoPoint2(posResGz)) { int succeedCount = 0; var cgResNew = cgRes.Clone(); for (int i = 0; i < confidenceCount; i++) { cgResNew.MainX = cgRes.MainX + RandomHelper.Normal(0, 200); cgResNew.MainY = cgRes.MainY + RandomHelper.Normal(0, 200); cgResNew.MainZ = cgRes.MainZ + RandomHelper.Normal(0, 200); cgResNew.Adja1X = cgRes.Adja1X + RandomHelper.Normal(0, 200); cgResNew.Adja1Y = cgRes.Adja1Y + RandomHelper.Normal(0, 200); cgResNew.Adja1Z = cgRes.Adja1Z + RandomHelper.Normal(0, 200); cgResNew.Dto1 = (cgRes.Dto1 * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; cgResNew.Dto2 = (cgRes.Dto2 * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; cgResNew.RefYbDto1 = (cgRes.RefYbDto1 * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; cgResNew.RefYbDto2 = (cgRes.RefYbDto2 * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; res = X2D1_GzPos(cgResNew, sRes); var posResNew = ConvertToGeoPoint(res); if (IsGeoPoint2(posResNew)) { var distance = PhysicsHelper.DistanceGeo((posResGz[0], posResGz[1], 0), (posResNew[0], posResNew[1], 0)); if (distance <= confidenceDistance) { succeedCount++; } } } posRes[6] = (int)(succeedCount / (double)confidenceCount * 100); } else { posRes[6] = 0; } } return posRes; } /// /// 两星一地无参,返回经度、纬度、高度、镜像点、置信度,数组长度为7 /// /// 参估结果 /// 站点信息 /// 是否计算置信度 /// public static double[] X2D1_PosNoRef(CgRes cgRes, StationRes sRes, bool CalcConfidence = false) { if (cgRes.Snr1 == 0 || cgRes.Snr2 == 0) { return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; } double[] mainSat = new double[3] { cgRes.MainX, cgRes.MainY, cgRes.MainZ }; double[] adjaSat = new double[3] { cgRes.Adja1X.Value, cgRes.Adja1Y.Value, cgRes.Adja1Z.Value }; double[] satStation = new double[3] { sRes.SatTxLon, sRes.SatTxLat, 0 }; double[] cdbStation = new double[3] { sRes.CdbTxLon.Value, sRes.CdbTxLat.Value, 0 }; double dto1 = cgRes.Dto1 / 1e6; double dtoCdb = cgRes.Dto2.Value / 1e6; double[] zone = new double[] { -85, 85, -180, 180 }; //定位区域 double[] res = new double[6]; X2D1_PosNoRef20240305_Core(mainSat, adjaSat, cdbStation, satStation, satStation, satStation, zone, dto1, dtoCdb, res); var posRes = ConvertToGeoPoint(res); if (CalcConfidence && IsGeoPoint2(posRes)) { int succeedCount = 0; for (int i = 0; i < confidenceCount; i++) { var mainSatNew = new double[3] { cgRes.MainX+RandomHelper.Normal(0,200), cgRes.MainY+RandomHelper.Normal(0,200), cgRes.MainZ+RandomHelper.Normal(0,200) }; var adjaSatNew = new double[3] { cgRes.Adja1X.Value+RandomHelper.Normal(0,200), cgRes.Adja1Y.Value+RandomHelper.Normal(0,200), cgRes.Adja1Z.Value+RandomHelper.Normal(0,200) }; var dto1New = (dto1 * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; var dtoCdbNew = (dtoCdb * 1e6 + RandomHelper.Normal(0, 1)) / 1e6; X2D1_PosNoRef20240305_Core(mainSatNew, adjaSatNew, cdbStation, satStation, satStation, satStation, zone, dto1New, dtoCdbNew, res); var posResNew = ConvertToGeoPoint(res); if (IsGeoPoint2(posResNew)) { var distance = PhysicsHelper.DistanceGeo((posRes[0], posRes[1], 0), (posResNew[0], posResNew[1], 0)); if (distance <= confidenceDistance) { succeedCount++; } } } posRes[6] = (int)(succeedCount / (double)confidenceCount * 100); } return posRes; } public static double[] X2D1_PosNoRef_ZL(CgRes cgRes, StationRes sRes, bool CalcConfidence = false) { if (cgRes.Snr1 == 0 || cgRes.Snr2 == 0) { return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; } var cdbStation = new double[2] { sRes.CdbTxLon.Value, sRes.CdbTxLat.Value }; var satStation = new double[2] { sRes.SatTxLon, sRes.SatTxLat }; var mainXYZ = new double[3] { cgRes.MainX, cgRes.MainY, cgRes.MainZ }; var adjaXYZ = new double[3] { cgRes.Adja1X.Value, cgRes.Adja1Y.Value, cgRes.Adja1Z.Value }; var dtoSxus = cgRes.Dto1; var dtoXdus = cgRes.Dto2.Value; var posRes = X2D1_PosNoRefCore_Zl(cdbStation, satStation, mainXYZ, adjaXYZ, dtoSxus, dtoXdus); posRes = ConvertToGeoPoint(posRes); if (CalcConfidence && IsGeoPoint2(posRes)) { int succeedCount = 0; int totalCount = 100; for (int i = 0; i < totalCount; i++) { var mainSatNew = new double[3] { cgRes.MainX+RandomHelper.Normal(0,500), cgRes.MainY+RandomHelper.Normal(0,500), cgRes.MainZ+RandomHelper.Normal(0,500) }; var adjaSatNew = new double[3] { cgRes.Adja1X.Value+RandomHelper.Normal(0,500), cgRes.Adja1Y.Value+RandomHelper.Normal(0,500), cgRes.Adja1Z.Value+RandomHelper.Normal(0,500) }; var dtoSxusNew = dtoSxus + RandomHelper.Normal(0, 1); var dtoXdusNew = dtoXdus + RandomHelper.Normal(0, 1); var posResNew = X2D1_PosNoRefCore_Zl(cdbStation, satStation, mainSatNew, adjaSatNew, dtoSxusNew, dtoXdusNew); posResNew = ConvertToGeoPoint(posResNew); if (IsGeoPoint2(posResNew)) { var distance = PhysicsHelper.DistanceGeo((posRes[0], posRes[1], 0), (posResNew[0], posResNew[1], 0)); if (distance <= confidenceDistance) { succeedCount++; } } } var confidence = (int)(succeedCount / (double)totalCount * 100); if (confidence == 0) confidence = 1; if (confidence == 100) confidence = 99; posRes[6] = confidence; } return posRes; } private static double[] X2D1_PosNoRefCore_Zl(double[] cdbStation, double[] satStation, double[] mainXYZ, double[] adjaXYZ, double dtoSxus, double dtoXdus) { double startLon = (int)cdbStation[0] - 10; double endLon = (int)cdbStation[0] + 10; double startLat = (int)cdbStation[1] - 10; double endLat = (int)cdbStation[1] + 10; List<(double, double, double, double)> list = new List<(double, double, double, double)>(); var recXYZ = PhysicsHelper.GeoToEcef((satStation[0], satStation[1], 0)); var cdbXYZ = PhysicsHelper.GeoToEcef((cdbStation[0], cdbStation[1], 0)); for (double lon = startLon; lon < endLon; lon += 0.0001) { int flag = 0; for (double lat = startLat; lat < endLat; lat += 0.0001) { var posXYZ = PhysicsHelper.GeoToEcef((lon, lat, 0)); //目标-主星-接收站的时间(us) var dt1 = PhysicsHelper.Dto(posXYZ, (mainXYZ[0], mainXYZ[1], mainXYZ[2]), recXYZ) * 1e6; //目标-超短站的时间(us) var dt3 = PhysicsHelper.Dto(posXYZ, cdbXYZ) * 1e6; var dto2 = Math.Abs(dt1 - dt3 - dtoXdus); if (dto2 > 400) { lat += 1; if (flag < 1) flag = 1; continue; } else if (dto2 > 200) { lat += 0.5; if (flag < 2) flag = 2; continue; } else if (dto2 > 100) { lat += 0.2; if (flag < 3) flag = 3; continue; } else if (dto2 > 10) { lat += 0.01; if (flag < 4) flag = 4; continue; } //else if (dto2 > 2) //{ // lat += 0.001; // if (flag < 5) // flag = 5; // continue; //} //目标-邻星-接收站的时间(us) var dt2 = PhysicsHelper.Dto(posXYZ, (adjaXYZ[0], adjaXYZ[1], adjaXYZ[2]), recXYZ) * 1e6; var dto1 = Math.Abs(dt1 - dt2 - dtoSxus); if (dto1 > 400) { lat += 1; if (flag < 1) flag = 1; continue; } if (dto1 > 200) { lat += 0.5; if (flag < 2) flag = 2; continue; } else if (dto1 > 100) { lat += 0.2; if (flag < 3) flag = 3; continue; } else if (dto1 > 10) { lat += 0.01; if (flag < 4) flag = 4; continue; } //else if (dto1 > 2) //{ // lat += 0.001; // if (flag < 5) // flag = 5; // continue; //} list.Add((lon, lat, dto1, dto2)); } if (flag == 1) lon += 1; else if (flag == 2) lon += 0.5; else if (flag == 3) lon += 0.2; else if (flag == 4) lon += 0.01; else if (flag == 5) lon += 0.001; } double[] posRes; var p1 = list.OrderBy(p => p.Item4 * p.Item4 + p.Item3 * p.Item3).FirstOrDefault(); if (p1 == default) { posRes = new double[6] { 999, 999, 0, 999, 999, 0 }; return posRes; } var p2 = list.Where(p => PhysicsHelper.DistanceGeo((p1.Item1, p1.Item2, 0), (p.Item1, p.Item2, 0)) > 10000).OrderBy(p => p.Item4 * p.Item4 + p.Item3 * p.Item3).FirstOrDefault(); if (p2 == default) posRes = new double[6] { p1.Item1, p1.Item2, 0, 999, 999, 0 }; else posRes = new double[6] { p1.Item1, p1.Item2, 0, p2.Item1, p2.Item2, 0 }; if (posRes[3] != 999) { var dis1 = PhysicsHelper.DistanceArcGeo((posRes[3], posRes[4]), (cdbStation[0], cdbStation[1])); var dis2 = PhysicsHelper.DistanceArcGeo((posRes[0], posRes[1]), (cdbStation[0], cdbStation[1])); if (dis1 < dis2) { var tmp1 = posRes[3]; var tmp2 = posRes[4]; posRes[3] = posRes[0]; posRes[4] = posRes[1]; posRes[0] = tmp1; posRes[1] = tmp2; } } return posRes; } private static double[] ConvertToGeoPoint(double[] res) { if (res == null || res.Length == 0) { return new double[7] { 999, 999, 0, 999, 999, 0, -1 }; } else if (res.Length != 3 && res.Length != 6) { throw new Exception("定位结果解析异常,长度不是3或6"); } var list = res.Select(p => Math.Round(p, 4)).ToArray(); if (list.Length == 3) { list = list.Concat(new double[4] { 999, 999, 0, -1 }).ToArray(); } list[2] = list[5] = 0;//高度 if (double.IsNaN(list[0]) || double.IsNaN(list[1])) { list[0] = list[1] = 999; } if (double.IsNaN(list[3]) || double.IsNaN(list[4])) { list[3] = list[4] = 999; } if (list[0] < -180 || list[0] > 180) { list[0] = list[1] = 999; } else if (list[1] < -90 || list[1] > 90) { list[0] = list[1] = 999; } if (list[3] < -180 || list[3] > 180) { list[3] = list[4] = 999; } else if (list[4] < -90 || list[4] > 90) { list[3] = list[4] = 999; } if (list.Length == 6) { var listNew = list.ToList(); listNew.Add(-1); list = listNew.ToArray(); } return list; } private static bool IsGeoPoint(double[] res) { if (res.Length != 3) return false; if (res[0] < -180 || res[0] > 180) { return false; } else if (res[1] < -90 || res[1] > 90) { return false; } return true; } private static bool IsGeoPoint2(double[] res) { if (res.Length < 3) return false; if (res[0] < -180 || res[0] > 180) { return false; } else if (res[1] < -90 || res[1] > 90) { return false; } return true; } } }