## 1 Introduction  
Time-series InSAR technology, which continuously captures phase information of surface deformation, has become an essential tool for urban subsidence monitoring and infrastructure health assessment. Its millimeter-level monitoring accuracy provides irreplaceable technical support for urban safety and maintenance [1, 2, 3, 4, 5]. Particularly for the health assessment of urban lifeline projects, high-efficiency, high-density, and high-precision time-series InSAR processing has become a prerequisite for accurate structural displacement inversion [6; 7; 8; 9]. Among various time-series InSAR methods, Permanent Scatterer InSAR (PSInSAR) stands out due to its unique phase-unwrapping strategy [10; 11; 12]. By selecting stable ground scatterers (PS points), it constructs a joint inversion model based on arc segments for the deformation rate and elevation residual estimation, effectively suppressing spatiotemporal decorrelation and atmospheric phase disturbances. This enables the retrieval of high-precision deformation and position information, making it particularly suitable for urban areas with dense artificial structures [13; 14]. In recent years, advancements in SAR satellite technology have led to three significant improvements: wide-swath observation capability (e.g., Sentinel-1 IW mode with a swath width of 250 km), high-frequency revisit cycles (e.g., LuTan-1 dual-satellite constellation achieving a 4-day revisit period [15]), and sub-meter spatial resolution (SAR spotlight mode resolution better than 1 m [16]). These improvements have drastically increased the density of PS points in urban areas, making it possible to resolve millimeter-scale deformations at individual building scales [17; 18; 19; 20]. However, they also pose significant challenges to PSInSAR processing. For example, in this study's case of Fuzhou's main urban area, a single TerraSAR-X stripmap image can contain up to 90,000 PS points per km\({}^{2}\), with the entire area comprising approximately 1.4 \(\times\) 10\({}^{7}\) PS points, forming a massive interferometric network. The deformation inversion of such a large-scale PS network requires solving a coefficient matrix with dimensions exceeding 10\({}^{7}\), far exceeding the capacity of a single computing unit. Traditional global PSInSAR methods, constrained by their O(N\({}^{3}\)) computational complexity and O(N\({}^{2}\)) memory requirements, struggle to efficiently process this enormous dataset. To address the challenge of large-scale, high-density PS point processing, current research has primarily focused on two dimensionality reduction strategies. The first approach is the hierarchical network method [21; 22; 23; 24; 25; 26], which draws inspiration from geodetic control networks with a \"coarse-to-fine\" framework. It first selects high-quality PS points using strict thresholds or sparse sampling to form a primary network. The remaining PS or Distributed Scatterer (DS) points are then used to densify the network into a secondary refined network. This method requires the primary network to maintain a high quality to ensure global network stability and connectivity [27]. To address potential ill-conditioning in the primary network, some studies [23; 24; 28] have employed ridge estimation techniques for stabilization. Meanwhile, other studies [29; 30] have introduced adaptive network reconstruction to link disconnected subnetworks into the largest connected component, thereby ensuring global connectivity. However, despite reducing the computational complexity, this approach still requires global storage of all the PS points and arc segment information, leaving memory overflow issues unresolved. The second approach is the partition-based processing method [31; 32; 33; 34], which divides the global processing task into smaller independent blocks, significantly reducing the memory usage and computational load per block. However, a major challenge is block merging.