无源侦察测向定位模型

#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <cmath>

// 假设数据结构如下
struct Measurement {
    double lambda_i; // 测向站经度(弧度)
    double phi_i;    // 测向站纬度(弧度)
    double theta_i;  // 测向方位角(弧度)，正北为0，顺时针方向
};

int main() {
    // 示例输入数据（需根据实际情况填写）
    // 假设有 n 次测量，每次给定测向站位置 (lambda_i, phi_i) 和方位 theta_i
    std::vector<Measurement> measurements = {
        // 示例数据，请根据实际情况修改
        { 1.234, 0.456, 0.523 },
        { 1.235, 0.457, 1.047 },
        { 1.237, 0.460, 1.570 }
    };

    int n = static_cast<int>(measurements.size());
    if (n < 2) {
        std::cerr << "至少需要两组测向数据才能求解！" << std::endl;
        return -1;
    }

    // 构建线性方程组 A * X = C
    // 对应关系：
    // A的每行: [cos(phi_i)*cos(theta_i), -sin(theta_i)]
    // C的每行: lambda_i * cos(phi_i)*cos(theta_i) - phi_i*sin(theta_i)
    //
    // X = [lambda_e, phi_e]^T 是未知量。

    Eigen::MatrixXd A(n, 2);
    Eigen::VectorXd C(n);

    for (int i = 0; i < n; ++i) {
        double lambda_i = measurements[i].lambda_i;
        double phi_i    = measurements[i].phi_i;
        double theta_i  = measurements[i].theta_i;

        double a_i = std::cos(phi_i)*std::cos(theta_i);
        double b_i = -std::sin(theta_i);
        double c_i = lambda_i * std::cos(phi_i)*std::cos(theta_i) - phi_i * std::sin(theta_i);

        A(i, 0) = a_i;
        A(i, 1) = b_i;
        C(i)    = c_i;
    }

    // 当 n=2 时可以直接 A.inverse()*C 求解
    // 当 n>2 时，需要使用最小二乘求解，即 (A^T A)X = A^T C
    // 使用Eigen提供的最小二乘求解器
    Eigen::VectorXd X;

    if (n == 2) {
        // 2x2可直接求逆
        Eigen::Matrix2d A2;
        Eigen::Vector2d C2;
        A2 << A(0,0), A(0,1),
              A(1,0), A(1,1);
        C2 << C(0), C(1);
        X = A2.inverse() * C2;
    } else {
        // n > 2，用最小二乘法
        // 正规方程: X = (A^T A)^{-1} A^T C
        // 或直接使用Eigen的colPivHouseholderQr等分解进行求解
        X = (A.transpose() * A).ldlt().solve(A.transpose() * C);
    }

    double lambda_e = X(0); // 求解出的辐射源经度
    double phi_e    = X(1); // 求解出的辐射源纬度

    // 输出结果
    std::cout << "求解得到的辐射源坐标：" << std::endl;
    std::cout << "lambda_e = " << lambda_e << " (弧度)" << std::endl;
    std::cout << "phi_e    = " << phi_e    << " (弧度)" << std::endl;

    // 若需要，可将弧度转换为度
    double lambda_e_deg = lambda_e * 180.0 / M_PI;
    double phi_e_deg = phi_e * 180.0 / M_PI;
    std::cout << "lambda_e (度) = " << lambda_e_deg << std::endl;
    std::cout << "phi_e    (度) = " << phi_e_deg  << std::endl;

    return 0;
}

#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <cmath>

// 测量数据结构
struct Measurement {
    double lambda_i; // 测向站经度 (弧度)
    double phi_i;    // 测向站纬度 (弧度)
    double theta_i;  // 测向方位角 (弧度)
};

int main() {
    // 示例测量数据，需根据实际情况修改
    std::vector<Measurement> measurements = {
        { 1.234, 0.456, 0.523 },
        { 1.235, 0.457, 1.047 },
        { 1.237, 0.460, 1.570 }
    };

    int n = static_cast<int>(measurements.size());
    if (n < 2) {
        std::cerr << "至少需要两组测向数据才能求解！" << std::endl;
        return -1;
    }

    Eigen::MatrixXd A(n, 2);
    Eigen::VectorXd C(n);

    // 构建方程组 A * X = C
    // A的每行对应: [cos(phi_i)*cos(theta_i), -sin(theta_i)]
    // C的每行对应: lambda_i*cos(phi_i)*cos(theta_i) - phi_i*sin(theta_i)
    for (int i = 0; i < n; ++i) {
        double lambda_i = measurements[i].lambda_i;
        double phi_i    = measurements[i].phi_i;
        double theta_i  = measurements[i].theta_i;

        double a_i = std::cos(phi_i)*std::cos(theta_i);
        double b_i = -std::sin(theta_i);
        double c_i = lambda_i * std::cos(phi_i)*std::cos(theta_i) - phi_i * std::sin(theta_i);

        A(i, 0) = a_i;
        A(i, 1) = b_i;
        C(i)    = c_i;
    }

    // 使用更加数值稳定的分解法求解最小二乘问题
    Eigen::VectorXd X = A.colPivHouseholderQr().solve(C);

    double lambda_e = X(0); // 辐射源经度 (弧度)
    double phi_e    = X(1); // 辐射源纬度 (弧度)

    std::cout << "求解得到的辐射源坐标(弧度)：" << std::endl;
    std::cout << "lambda_e = " << lambda_e << std::endl;
    std::cout << "phi_e    = " << phi_e << std::endl;

    // 若需要度制输出
    double lambda_e_deg = lambda_e * 180.0 / M_PI;
    double phi_e_deg = phi_e * 180.0 / M_PI;

    std::cout << "辐射源坐标(度)：" << std::endl;
    std::cout << "lambda_e (deg) = " << lambda_e_deg << std::endl;
    std::cout << "phi_e    (deg) = " << phi_e_deg << std::endl;

    return 0;
}

import React, { useState } from 'react';
import { Card, CardHeader, CardTitle, CardContent } from '@/components/ui/card';
import { Button } from '@/components/ui/button';
import { Input } from '@/components/ui/input';

const RadiationCalculator = () => {
  const [measurements, setMeasurements] = useState([
    { lambda_i: '', phi_i: '', theta_i: '' }
  ]);
  const [actualSource, setActualSource] = useState({
    lambda_e: '',
    phi_e: ''
  });
  const [result, setResult] = useState(null);

  // 计算两点间的距离（米）
  const calculateDistance = (lat1, lon1, lat2, lon2) => {
    const R = 6371000; // 地球平均半径（米）
    const phi1 = lat1 * Math.PI / 180;
    const phi2 = lat2 * Math.PI / 180;
    const deltaPhi = (lat2 - lat1) * Math.PI / 180;
    const deltaLambda = (lon2 - lon1) * Math.PI / 180;

    const a = Math.sin(deltaPhi/2) * Math.sin(deltaPhi/2) +
              Math.cos(phi1) * Math.cos(phi2) *
              Math.sin(deltaLambda/2) * Math.sin(deltaLambda/2);
    const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    const distance = R * c;

    return distance;
  };

  // 计算所有站点之间的距离
  const calculateAllDistances = (sourcePos, stations) => {
    const distances = {
      sourceToStations: [], // 辐射源到各测向站的距离
      betweenStations: [] // 测向站之间的距离
    };

    // 计算辐射源到各测向站的距离
    stations.forEach((station, idx) => {
      const distance = calculateDistance(
        sourcePos.phi,
        sourcePos.lambda,
        parseFloat(station.phi_i),
        parseFloat(station.lambda_i)
      );
      distances.sourceToStations.push({
        stationIndex: idx + 1,
        distance: distance / 1000 // 转换为公里
      });
    });

    // 计算测向站之间的距离
    for (let i = 0; i < stations.length; i++) {
      for (let j = i + 1; j < stations.length; j++) {
        const distance = calculateDistance(
          parseFloat(stations[i].phi_i),
          parseFloat(stations[i].lambda_i),
          parseFloat(stations[j].phi_i),
          parseFloat(stations[j].lambda_i)
        );
        distances.betweenStations.push({
          station1: i + 1,
          station2: j + 1,
          distance: distance / 1000 // 转换为公里
        });
      }
    }

    return distances;
  };

  const addMeasurement = () => {
    setMeasurements([...measurements, { lambda_i: '', phi_i: '', theta_i: '' }]);
  };

  const removeMeasurement = (index) => {
    if (measurements.length > 1) {
      const newMeasurements = measurements.filter((_, i) => i !== index);
      setMeasurements(newMeasurements);
    }
  };

  const handleMeasurementChange = (index, field, value) => {
    const newMeasurements = [...measurements];
    newMeasurements[index][field] = value;
    setMeasurements(newMeasurements);
  };

  const calculatePosition = () => {
    // 将度数转换为弧度
    const measurementsRad = measurements.map(m => ({
      lambda_i: parseFloat(m.lambda_i) * Math.PI / 180,
      phi_i: parseFloat(m.phi_i) * Math.PI / 180,
      theta_i: parseFloat(m.theta_i) * Math.PI / 180
    }));

    const n = measurementsRad.length;
    if (n < 2) {
      alert('至少需要两组测向数据！');
      return;
    }

    // 构建方程组
    const A = new Array(n).fill(0).map(() => new Array(2).fill(0));
    const C = new Array(n).fill(0);

    for (let i = 0; i < n; i++) {
      const { lambda_i, phi_i, theta_i } = measurementsRad[i];
      const cos_phi = Math.cos(phi_i);
      const cos_theta = Math.cos(theta_i);
      const sin_theta = Math.sin(theta_i);

      A[i][0] = cos_phi * cos_theta;
      A[i][1] = -sin_theta;
      C[i] = lambda_i * cos_phi * cos_theta - phi_i * sin_theta;
    }

    // 使用最小二乘法求解
    const AT = A[0].map((_, colIndex) => A.map(row => row[colIndex]));
    const ATA = AT.map(row => {
      return A[0].map((_, j) => {
        return row.reduce((sum, _, k) => sum + row[k] * A[k][j], 0);
      });
    });
    const ATC = AT.map(row => {
      return row.reduce((sum, _, i) => sum + row[i] * C[i], 0);
    });

    const det = ATA[0][0] * ATA[1][1] - ATA[0][1] * ATA[1][0];
    const lambda_e = (ATA[1][1] * ATC[0] - ATA[0][1] * ATC[1]) / det;
    const phi_e = (-ATA[1][0] * ATC[0] + ATA[0][0] * ATC[1]) / det;

    // 转换回度数
    const lambda_e_deg = lambda_e * 180 / Math.PI;
    const phi_e_deg = phi_e * 180 / Math.PI;

    // 计算与实际位置的偏差
    const actual_lambda = parseFloat(actualSource.lambda_e);
    const actual_phi = parseFloat(actualSource.phi_e);

    const lambda_error = Math.abs(lambda_e_deg - actual_lambda);
    const phi_error = Math.abs(phi_e_deg - actual_phi);

    // 计算距离偏差
    const distance_error = calculateDistance(
      actual_phi,
      actual_lambda,
      phi_e_deg,
      lambda_e_deg
    );

    // 计算所有站点之间的距离
    const calculatedDistances = calculateAllDistances(
      { lambda: lambda_e_deg, phi: phi_e_deg },
      measurements
    );
    
    const actualDistances = calculateAllDistances(
      { lambda: actual_lambda, phi: actual_phi },
      measurements
    );

    setResult({
      calculated: {
        lambda: lambda_e_deg.toFixed(4),
        phi: phi_e_deg.toFixed(4)
      },
      error: {
        lambda: lambda_error.toFixed(4),
        phi: phi_error.toFixed(4),
        distance: (distance_error / 1000).toFixed(3)  // 转换为公里
      },
      distances: {
        calculated: calculatedDistances,
        actual: actualDistances
      }
    });
  };

  return (
    <Card className="w-full max-w-4xl">
      <CardHeader>
        <CardTitle>辐射源定位计算器</CardTitle>
      </CardHeader>
      <CardContent className="space-y-6">
        <div className="space-y-4">
          <div className="font-medium">实际辐射源位置（度）：</div>
          <div className="grid grid-cols-2 gap-4">
            <div>
              <label className="block text-sm mb-1">辐射源经度 (λe)：</label>
              <Input
                type="number"
                value={actualSource.lambda_e}
                onChange={(e) => setActualSource({...actualSource, lambda_e: e.target.value})}
                placeholder="输入实际经度"
              />
            </div>
            <div>
              <label className="block text-sm mb-1">辐射源纬度 (φe)：</label>
              <Input
                type="number"
                value={actualSource.phi_e}
                onChange={(e) => setActualSource({...actualSource, phi_e: e.target.value})}
                placeholder="输入实际纬度"
              />
            </div>
          </div>
        </div>

        <div className="space-y-4">
          <div className="font-medium">测向站数据（度）：</div>
          {measurements.map((m, index) => (
            <div key={index} className="grid grid-cols-3 gap-4 items-end">
              <div>
                <label className="block text-sm mb-1">测向站经度 (λi)：</label>
                <Input
                  type="number"
                  value={m.lambda_i}
                  onChange={(e) => handleMeasurementChange(index, 'lambda_i', e.target.value)}
                  placeholder="输入测向站经度"
                />
              </div>
              <div>
                <label className="block text-sm mb-1">测向站纬度 (φi)：</label>
                <Input
                  type="number"
                  value={m.phi_i}
                  onChange={(e) => handleMeasurementChange(index, 'phi_i', e.target.value)}
                  placeholder="输入测向站纬度"
                />
              </div>
              <div>
                <label className="block text-sm mb-1">方位角 (θi)：</label>
                <div className="flex gap-2">
                  <Input
                    type="number"
                    value={m.theta_i}
                    onChange={(e) => handleMeasurementChange(index, 'theta_i', e.target.value)}
                    placeholder="输入方位角"
                  />
                  {measurements.length > 1 && (
                    <Button 
                      variant="destructive"
                      className="px-3"
                      onClick={() => removeMeasurement(index)}
                    >
                      删除
                    </Button>
                  )}
                </div>
              </div>
            </div>
          ))}
        </div>

        <div className="space-x-4">
          <Button onClick={addMeasurement}>添加测向站</Button>
          <Button onClick={calculatePosition}>计算位置</Button>
        </div>

        {result && (
          <div className="mt-6 p-4 bg-gray-100 rounded-lg">
            <h3 className="font-medium mb-2">计算结果：</h3>
            <div className="space-y-4">
              <div>
                <p className="font-medium">计算得到的位置：</p>
                <p>经度 (λe) = {result.calculated.lambda}°</p>
                <p>纬度 (φe) = {result.calculated.phi}°</p>
              </div>
              
              <div>
                <p className="font-medium">与实际位置的偏差：</p>
                <p>经度偏差 = {result.error.lambda}°</p>
                <p>纬度偏差 = {result.error.phi}°</p>
                <p>距离偏差 = {result.error.distance} 公里</p>
              </div>

              <div>
                <p className="font-medium">计算位置到各测向站的距离：</p>
                {result.distances.calculated.sourceToStations.map((dist, idx) => (
                  <p key={idx}>
                    到测向站{dist.stationIndex}: {dist.distance.toFixed(3)} 公里
                  </p>
                ))}
              </div>

              <div>
                <p className="font-medium">测向站之间的距离：</p>
                {result.distances.calculated.betweenStations.map((dist, idx) => (
                  <p key={idx}>
                    测向站{dist.station1}到测向站{dist.station2}: {dist.distance.toFixed(3)} 公里
                  </p>
                ))}
              </div>
            </div>
          </div>
        )}
      </CardContent>
    </Card>
  );
};

export default RadiationCalculator;